VesuvioCorrections.py

# Mantid Repository : https://github.com/mantidproject/mantid
#
# Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
#   NScD Oak Ridge National Laboratory, European Spallation Source,
#   Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
# SPDX - License - Identifier: GPL - 3.0 +
# pylint: disable=no-init, too-many-instance-attributes
from mantid.kernel import *
from mantid.api import *
from vesuvio.base import VesuvioBase, TableWorkspaceDictionaryFacade
from vesuvio.fitting import parse_fit_options
from vesuvio.instrument import VESUVIO
import mantid.simpleapi as ms

import math


# ----------------------------------------------------------------------------------------


def create_cuboid_xml(height, width, depth):
    """
    Create the XML string to describe a cuboid of the given dimensions

    @param height Height in metres (Y coordinate)
    @param width Width in metres (X coordinate)
    @param depth Depth in metres (Z coordinate)
    """
    half_height, half_width, half_thick = 0.5 * height, 0.5 * width, 0.5 * depth
    xml_str = \
        " <cuboid id=\"sample-shape\"> " \
        + "<left-front-bottom-point " \
        + "x=\"%f\" y=\"%f\" z=\"%f\" /> " % (half_width, -half_height, half_thick) \
        + "<left-front-top-point " \
        + "x=\"%f\" y=\"%f\" z=\"%f\" /> " % (half_width, half_height, half_thick) \
        + "<left-back-bottom-point " \
        + "x=\"%f\" y=\"%f\" z=\"%f\" /> " % (half_width, -half_height, -half_thick) \
        + "<right-front-bottom-point " \
        + "x=\"%f\" y=\"%f\" z=\"%f\" /> " % (-half_width, -half_height, half_thick) \
        + "</cuboid>"
    return xml_str


# ----------------------------------------------------------------------------------------


class VesuvioCorrections(VesuvioBase):
    _input_ws = None
    _output_ws = None
    _correction_workspaces = None
    _linear_fit_table = None
    _correction_wsg = None
    _corrected_wsg = None
    _container_ws = None
    _spec_idx = None

    # ------------------------------------------------------------------------------

    def summary(self):
        return "Apply post fitting steps to vesuvio data"

    def category(self):
        return "Inelastic\\Indirect\\Vesuvio"

    # ------------------------------------------------------------------------------

    # pylint: disable=too-many-locals
    def PyInit(self):

        # -------------------------------------------------------------------------------------------

        # Input Property setup

        self.declareProperty(MatrixWorkspaceProperty("InputWorkspace", "",
                                                     direction=Direction.Input),
                             doc="Input TOF workspace")

        self.declareProperty("WorkspaceIndex", 0,
                             doc="Index of spectrum to calculate corrections for")

        self.declareProperty(ITableWorkspaceProperty("FitParameters", "",
                                                     direction=Direction.Input,
                                                     optional=PropertyMode.Optional),
                             doc="Table containing the calculated fit parameters"
                                 "for the data in the workspace")

        input_group = "Input Options"
        self.setPropertyGroup("InputWorkspace", input_group)
        self.setPropertyGroup("WorkspaceIndex", input_group)
        self.setPropertyGroup("FitParameters", input_group)

        # -------------------------------------------------------------------------------------------

        # Mass Property setup

        float_length_validator = FloatArrayLengthValidator()
        float_length_validator.setLengthMin(1)
        self.declareProperty(FloatArrayProperty("Masses", float_length_validator),
                             doc="Mass values for fitting")

        self.declareProperty(PropertyManagerProperty("MassIndexToSymbolMap", {},
                                                     direction=Direction.Input),
                             doc="A map from the index of the mass in the Masses"
                                 " property to a chemical symbol.")

        self.declareProperty("MassProfiles", "", StringMandatoryValidator(),
                             doc="Functions used to approximate mass profile. "
                                 "The format is "
                                 "function=Function1Name,param1=val1,param2=val2;"
                                 "function=Function2Name,param3=val3,param4=val4")

        mass_group = "Mass Options"
        self.setPropertyGroup("Masses", mass_group)
        self.setPropertyGroup("MassIndexToSymbolMap", mass_group)
        self.setPropertyGroup("MassProfiles", mass_group)

        # -------------------------------------------------------------------------------------------

        # Mass Constraints Property setup

        self.declareProperty("IntensityConstraints", "",
                             doc="A semi-colon separated list of intensity "
                                 "constraints defined as lists e.g "
                                 "[0,1,0,-4];[1,0,-2,0]")

        self.declareProperty(PropertyManagerProperty("HydrogenConstraints", {},
                                                     direction=Direction.Input),
                             doc="Constraints used to approximate the intensity of"
                                 " the hydrogen peak in back-scattering spectra for"
                                 " multiple scattering corrections.")

        mass_constraints_group = "Mass Constraints"
        self.setPropertyGroup("IntensityConstraints", mass_constraints_group)
        self.setPropertyGroup("HydrogenConstraints", mass_constraints_group)

        # -------------------------------------------------------------------------------------------

        # Container Property setup

        self.declareProperty(MatrixWorkspaceProperty("ContainerWorkspace", "",
                                                     direction=Direction.Input,
                                                     optional=PropertyMode.Optional),
                             doc="Container workspace in TOF")

        self.declareProperty("ContainerScale", 0.0,
                             doc="Scale factor to apply to container, set to 0 for "
                                 "automatic scale based on linear fit")

        container_group = "Container Options"
        self.setPropertyGroup("ContainerWorkspace", container_group)
        self.setPropertyGroup("ContainerScale", container_group)

        # -------------------------------------------------------------------------------------------

        # Gamma Background Property setup

        self.declareProperty("GammaBackground", True, direction=Direction.Input,
                             doc="If true, correct for the gamma background")

        self.declareProperty("GammaBackgroundScale", 0.0,
                             doc="Scale factor to apply to gamma background, set to 0 "
                                 "for automatic scale based on linear fit")

        gamma_group = "Gamma Correction Options"
        self.setPropertyGroup("GammaBackground", gamma_group)
        self.setPropertyGroup("GammaBackgroundScale", gamma_group)

        gamma_enabled = VisibleWhenProperty('GammaBackground', PropertyCriterion.IsEqualTo, "1")
        self.setPropertySettings("GammaBackgroundScale", gamma_enabled)

        # -------------------------------------------------------------------------------------------

        # Multiple Scattering Property setup

        self.declareProperty("MultipleScattering", True, direction=Direction.Input,
                             doc="If true, correct for the effects of multiple scattering")

        self.declareProperty("BeamRadius", 2.5,
                             doc="Radius of beam in cm")

        self.declareProperty("SampleHeight", 5.0,
                             doc="Height of sample in cm")

        self.declareProperty("SampleWidth", 5.0,
                             doc="Width of sample in cm")

        self.declareProperty("SampleDepth", 5.0,
                             doc="Depth of sample in cm")

        self.declareProperty("SampleDensity", 1.0,
                             doc="Sample density in g/cm^3")

        self.declareProperty("Seed", 123456789,
                             doc="")

        self.declareProperty("NumScatters", 3,
                             doc="")

        self.declareProperty("NumRuns", 10,
                             doc="")

        self.declareProperty("NumEvents", 50000,
                             doc="Number of neutron events")

        self.declareProperty("SmoothNeighbours", 3,
                             doc="")

        ms_group = "Multiple Scattering Options"
        self.setPropertyGroup("MultipleScattering", ms_group)
        self.setPropertyGroup("BeamRadius", ms_group)
        self.setPropertyGroup("SampleHeight", ms_group)
        self.setPropertyGroup("SampleWidth", ms_group)
        self.setPropertyGroup("SampleDepth", ms_group)
        self.setPropertyGroup("SampleDensity", ms_group)
        self.setPropertyGroup("Seed", ms_group)
        self.setPropertyGroup("NumScatters", ms_group)
        self.setPropertyGroup("NumRuns", ms_group)
        self.setPropertyGroup("NumEvents", ms_group)
        self.setPropertyGroup("SmoothNeighbours", ms_group)

        ms_enabled = VisibleWhenProperty('MultipleScattering', PropertyCriterion.IsEqualTo, "1")
        self.setPropertySettings("BeamRadius", ms_enabled)
        self.setPropertySettings("SampleHeight", ms_enabled)
        self.setPropertySettings("SampleWidth", ms_enabled)
        self.setPropertySettings("SampleDepth", ms_enabled)
        self.setPropertySettings("SampleDensity", ms_enabled)
        self.setPropertySettings("Seed", ms_enabled)
        self.setPropertySettings("NumScatters", ms_enabled)
        self.setPropertySettings("NumRuns", ms_enabled)
        self.setPropertySettings("NumEvents", ms_enabled)
        self.setPropertySettings("SmoothNeighbours", ms_enabled)

        # Disable hydrogen constraints when there is no multiple scattering
        self.setPropertySettings("HydrogenConstraints", ms_enabled)

        # -------------------------------------------------------------------------------------------

        # Outputs Property setup

        self.declareProperty(WorkspaceGroupProperty("CorrectionWorkspaces", "",
                                                    direction=Direction.Output,
                                                    optional=PropertyMode.Optional),
                             doc="Workspace group containing correction intensities "
                                 "for each correction")

        self.declareProperty(WorkspaceGroupProperty("CorrectedWorkspaces", "",
                                                    direction=Direction.Output,
                                                    optional=PropertyMode.Optional),
                             doc="Workspace group containing individual corrections "
                                 "applied to raw data")

        self.declareProperty(ITableWorkspaceProperty("LinearFitResult", "",
                                                     direction=Direction.Output,
                                                     optional=PropertyMode.Optional),
                             doc="Table workspace containing the fit parameters used to"
                                 "linearly fit the corrections to the data")

        self.declareProperty(MatrixWorkspaceProperty("OutputWorkspace", "",
                                                     direction=Direction.Output),
                             doc="The name of the output workspace")

        output_group = "Output Options"
        self.setPropertyGroup("CorrectionWorkspaces", output_group)
        self.setPropertyGroup("CorrectedWorkspaces", output_group)
        self.setPropertyGroup("LinearFitResult", output_group)
        self.setPropertyGroup("OutputWorkspace", output_group)

    # ------------------------------------------------------------------------------

    def validateInputs(self):
        self._get_properties()
        errors = dict()

        if self.getProperty("FitParameters").value is None:
            errors["FitParameters"] = "Corrections require a " \
                                      + "set of parameters from a fit of the data"

        return errors

    # ------------------------------------------------------------------------------

    def _get_properties(self):
        self._input_ws = self.getProperty("InputWorkspace").value
        self._container_ws = self.getPropertyValue("ContainerWorkspace")
        self._spec_idx = self.getProperty("WorkspaceIndex").value
        self._output_ws = self.getPropertyValue("OutputWorkspace")
        self._correction_wsg = self.getPropertyValue("CorrectionWorkspaces")
        self._corrected_wsg = self.getPropertyValue("CorrectedWorkspaces")
        self._linear_fit_table = self.getPropertyValue("LinearFitResult")
        self._masses = self.getProperty("Masses").value
        self._index_to_symbol_map = self.getProperty("MassIndexToSymbolMap").value
        self._hydrogen_constraints = self.getProperty("HydrogenConstraints").value
        spec_no = self._input_ws.getSpectrum(self._spec_idx).getSpectrumNo()
        back_spectra = VESUVIO().backward_spectra
        self._back_scattering = back_spectra[0] <= spec_no <= back_spectra[1]

    # ------------------------------------------------------------------------------

    def PyExec(self):
        ms.ExtractSingleSpectrum(InputWorkspace=self._input_ws,
                                 OutputWorkspace=self._output_ws,
                                 WorkspaceIndex=self._spec_idx)

        # Performs corrections
        self._define_corrections()

        # The workspaces to fit for correction scale factors
        fit_corrections = [wks for wks in self._correction_workspaces if 'MultipleScattering' not in wks]

        # Perform fitting of corrections
        fixed_params = {}

        fixed_gamma_factor = self.getProperty("GammaBackgroundScale").value
        if fixed_gamma_factor != 0.0 and not self._back_scattering:
            fixed_params['GammaBackground'] = fixed_gamma_factor

        fixed_container_scale = self.getProperty("ContainerScale").value
        if fixed_container_scale != 0.0:
            fixed_params['Container'] = fixed_container_scale

        params_ws = self._fit_corrections(fit_corrections, self._linear_fit_table, **fixed_params)
        self.setProperty("LinearFitResult", params_ws)

        # Scale gamma background
        if self.getProperty("GammaBackground").value and not self._back_scattering:
            gamma_correct_ws = self._get_correction_workspace('GammaBackground')[1]
            gamma_factor = self._get_correction_scale_factor('GammaBackground',
                                                             fit_corrections, params_ws)
            ms.Scale(InputWorkspace=gamma_correct_ws,
                     OutputWorkspace=gamma_correct_ws,
                     Factor=gamma_factor)

        # Scale multiple scattering
        if self.getProperty("MultipleScattering").value:
            # Use factor of total scattering as this includes single and multiple scattering
            multi_scatter_correct_ws = self._get_correction_workspace('MultipleScattering')[1]
            total_scatter_correct_ws = self._get_correction_workspace('TotalScattering')[1]
            total_scatter_factor = self._get_correction_scale_factor('TotalScattering',
                                                                     fit_corrections, params_ws)
            ms.Scale(InputWorkspace=multi_scatter_correct_ws,
                     OutputWorkspace=multi_scatter_correct_ws,
                     Factor=total_scatter_factor)
            ms.Scale(InputWorkspace=total_scatter_correct_ws,
                     OutputWorkspace=total_scatter_correct_ws,
                     Factor=total_scatter_factor)

        # Scale by container
        if self._container_ws != "":
            container_correct_ws = self._get_correction_workspace('Container')[1]
            container_factor = self._get_correction_scale_factor('Container',
                                                                 fit_corrections, params_ws)
            ms.Scale(InputWorkspace=container_correct_ws,
                     OutputWorkspace=container_correct_ws,
                     Factor=container_factor)

        # Calculate and output corrected workspaces as a WorkspaceGroup
        if self._corrected_wsg != "":
            corrected_workspaces = [ws_name.replace(self._correction_wsg, self._corrected_wsg)
                                    for ws_name in self._correction_workspaces]
            for corrected, correction in zip(corrected_workspaces, self._correction_workspaces):
                ms.Minus(LHSWorkspace=self._output_ws,
                         RHSWorkspace=correction,
                         OutputWorkspace=corrected)
            ms.GroupWorkspaces(InputWorkspaces=corrected_workspaces,
                               OutputWorkspace=self._corrected_wsg)
            self.setProperty("CorrectedWorkspaces", self._corrected_wsg)

        # Apply corrections
        for correction in self._correction_workspaces:
            if 'TotalScattering' not in correction:
                ms.Minus(LHSWorkspace=self._output_ws,
                         RHSWorkspace=correction,
                         OutputWorkspace=self._output_ws)

        self.setProperty("OutputWorkspace", self._output_ws)

        # Remove correction workspaces if they are no longer required
        if self._correction_wsg == "":
            for wksp in self._correction_workspaces:
                ms.DeleteWorkspace(wksp)

                # ------------------------------------------------------------------------------

    def _define_corrections(self):
        """
        Defines all the corrections that are required
        """
        self._correction_workspaces = list()

        if self._container_ws != "":
            container_name = str(self._correction_wsg) + "_Container"
            self._container_ws = ms.ExtractSingleSpectrum(InputWorkspace=self._container_ws,
                                                          OutputWorkspace=container_name,
                                                          WorkspaceIndex=self._spec_idx)
            self._correction_workspaces.append(self._container_ws.name())

        # Do gamma correction
        if self.getProperty("GammaBackground").value and not self._back_scattering:
            self._correction_workspaces.append(self._gamma_correction())

        # Do multiple scattering correction
        if self.getProperty("MultipleScattering").value:
            self._correction_workspaces.extend(self._ms_correction())

        # Output correction workspaces as a WorkspaceGroup
        if self._correction_wsg != "":
            ms.GroupWorkspaces(InputWorkspaces=self._correction_workspaces,
                               OutputWorkspace=self._correction_wsg)
            self.setProperty("CorrectionWorkspaces", self._correction_wsg)

            # ------------------------------------------------------------------------------

    def _fit_corrections(self, fit_workspaces, param_table_name, **fixed_parameters):
        functions = []

        for idx, wsn in enumerate(fit_workspaces):
            tie = ''
            for param, value in fixed_parameters.items():
                if param in wsn:
                    tie = 'Scaling=%f,' % value
            function_str = "name=TabulatedFunction,Workspace=%s," % wsn \
                           + "ties=(%sShift=0,XScaling=1)," % tie \
                           + "constraints=(Scaling>=0.0)"
            functions.append(function_str)

            logger.notice('Corrections scale fit index %d is %s' % (idx, wsn))

        fit = AlgorithmManager.create("Fit")
        fit.initialize()
        fit.setChild(True)
        fit.setLogging(True)
        fit.setProperty("Function", ";".join(functions))
        fit.setProperty("InputWorkspace", self._output_ws)
        fit.setProperty("Output", param_table_name)
        fit.setProperty("CreateOutput", True)
        fit.execute()

        return fit.getProperty('OutputParameters').value

    # ------------------------------------------------------------------------------

    def _get_correction_workspace(self, correction_name, corrections=None):
        if corrections is None:
            corrections = self._correction_workspaces

        for idx, ws_name in enumerate(corrections):
            if correction_name in ws_name:
                return idx, ws_name

        return None, None

    # ------------------------------------------------------------------------------

    def _get_correction_scale_factor(self, correction_name, corrections, params_ws):
        index = self._get_correction_workspace(correction_name, corrections)[0]
        if index is None:
            raise RuntimeError('No workspace for given correction')

        params_dict = TableWorkspaceDictionaryFacade(params_ws)

        if len(corrections) > 1:
            scale_param_name = 'f%d.Scaling' % index
        else:
            scale_param_name = 'Scaling'

        return params_dict[scale_param_name]

    # ------------------------------------------------------------------------------

    def _gamma_correction(self):
        correction_background_ws = str(self._correction_wsg) + "_GammaBackground"

        fit_opts = parse_fit_options(mass_values=self._masses,
                                     profile_strs=self.getProperty("MassProfiles").value,
                                     constraints_str=self.getProperty("IntensityConstraints").value)
        params_dict = TableWorkspaceDictionaryFacade(self.getProperty("FitParameters").value)
        func_str = fit_opts.create_function_str(params_dict)

        ms.VesuvioCalculateGammaBackground(InputWorkspace=self._output_ws,
                                           ComptonFunction=func_str,
                                           BackgroundWorkspace=correction_background_ws,
                                           CorrectedWorkspace='__corrected_dummy')
        ms.DeleteWorkspace('__corrected_dummy')

        return correction_background_ws

    # ------------------------------------------------------------------------------

    def _ms_correction(self):
        """
        Calculates the contributions from multiple scattering
        on the input data from the set of given options
        """
        params_dict = TableWorkspaceDictionaryFacade(self.getProperty("FitParameters").value)

        atom_props = list()
        intensities = list()

        contains_hydrogen = False

        i = 0

        for idx, mass in enumerate(self._masses):

            if str(idx) in self._index_to_symbol_map:
                symbol = self._index_to_symbol_map[str(idx)].value
            else:
                symbol = None

            if symbol == 'H' and self._back_scattering:
                contains_hydrogen = True
                continue

            intensity_prop = 'f%d.Intensity' % i
            c0_prop = 'f%d.C_0' % i

            if intensity_prop in params_dict:
                intensity = params_dict[intensity_prop]
            elif c0_prop in params_dict:
                intensity = params_dict[c0_prop]
            else:
                i = i + 1
                continue

            # The program DINSMS_BATCH uses those sample parameters together with the sigma divided
            # by the sum absolute of scattering intensities for each detector (detector bank),
            # sigma/int_sum
            # Thus:
            # intensity = intensity/intensity_sum

            # In the thin sample limit, 1-exp(-n*dens*sigma) ~ n*dens*sigma, effectively the same
            # scattering power (ratio of double to single scatt.)  is obtained either by using
            # relative intensities ( sigma/int_sum ) or density divided by the total intensity
            # However, in the realistic case of thick sample, the SampleDensity, dens,  must be
            # obtained by iterative numerical solution of the Eq:
            # 1-exp(-n*dens*sigma) = measured scattering power of the sample.
            # For this, a program like THICK must be used.
            # The program THICK also uses sigma/int_sum to be consistent with the prgram
            # DINSMS_BATCH

            width_prop = 'f%d.Width' % i
            sigma_x_prop = 'f%d.SigmaX' % i
            sigma_y_prop = 'f%d.SigmaY' % i
            sigma_z_prop = 'f%d.SigmaZ' % i

            if width_prop in params_dict:
                width = params_dict['f%d.Width' % i]
            elif sigma_x_prop in params_dict:
                sigma_x = float(params_dict[sigma_x_prop])
                sigma_y = float(params_dict[sigma_y_prop])
                sigma_z = float(params_dict[sigma_z_prop])
                width = math.sqrt((sigma_x ** 2 + sigma_y ** 2 + sigma_z ** 2) / 3.0)
            else:
                i = i + 1
                continue

            atom_props.append(mass)
            atom_props.append(intensity)
            atom_props.append(width)
            intensities.append(intensity)

            # Check for NoneType is necessary as hydrogen constraints are
            # stored in a C++ PropertyManager object, not a dict; call to
            # __contains__ must match the C++ signature.
            if self._back_scattering and symbol is not None and symbol in self._hydrogen_constraints:
                self._hydrogen_constraints[symbol].value['intensity'] = intensity

            i = i + 1

        if self._back_scattering and contains_hydrogen:
            material_builder = MaterialBuilder()
            hydrogen = material_builder.setFormula('H').build()
            hydrogen_intensity = \
                self._calculate_hydrogen_intensity(hydrogen, self._hydrogen_constraints)
            hydrogen_width = 5
            atom_props.append(hydrogen.relativeMolecularMass())
            atom_props.append(hydrogen_intensity)
            atom_props.append(hydrogen_width)
            intensities.append(hydrogen_intensity)

        intensity_sum = sum(intensities)

        # Create the sample shape
        # Input dimensions are expected in CM
        ms.CreateSampleShape(InputWorkspace=self._output_ws,
                             ShapeXML=create_cuboid_xml(self.getProperty("SampleHeight").value / 100.,
                                                        self.getProperty("SampleWidth").value / 100.,
                                                        self.getProperty("SampleDepth").value / 100.))

        # Massage options into how algorithm expects them
        total_scatter_correction = str(self._correction_wsg) + "_TotalScattering"
        multi_scatter_correction = str(self._correction_wsg) + "_MultipleScattering"

        # Calculation
        # In the thin sample limit, 1-exp(-n*dens*sigma) ~ n*dens*sigma, effectively the same
        # scattering power(ratio of double to single scatt.)  is obtained either by using relative
        # intensities ( sigma/int_sum )or density divided by the total intensity.
        # However, in the realistic case of thick sample, the SampleDensity, dens,  must be
        # obtained by iterative numerical solution of the Eq:
        # 1-exp(-n*dens*sigma) = measured scattering power of the sample.
        # For this, a program like THICK must be used.
        # The program THICK also uses sigma/int_sum to be consistent with the prgram DINSMS_BATCH
        # The algorithm VesuvioCalculateMs called by the algorithm VesuvioCorrections takes the
        # parameter AtomicProperties with the absolute intensities, contraty to DINSMS_BATCH which
        # takes in relative intensities.
        # To compensate for this, the thickness parameter, dens (SampleDensity),  is divided in by
        # the sum of absolute intensities in VesuvioCorrections before being passed to
        # VesuvioCalculateMs.
        # Then, for the modified VesuvioCorrection algorithm one can use the thickenss parameter is
        # as is from the THICK command, i.e. 43.20552
        # This works, however, only in the thin sample limit, contrary to the THICK program. Thus,
        # for some detectors (detector banks) the SampleDensiy parameter may be over(under)
        # estimated.

        ms.VesuvioCalculateMS(InputWorkspace=self._output_ws,
                              NoOfMasses=int(len(atom_props) / 3),
                              SampleDensity=self.getProperty("SampleDensity").value / intensity_sum,
                              AtomicProperties=atom_props,
                              BeamRadius=self.getProperty("BeamRadius").value,
                              NumEventsPerRun=self.getProperty("NumEvents").value,
                              TotalScatteringWS=total_scatter_correction,
                              MultipleScatteringWS=multi_scatter_correction)

        # Smooth the output
        smooth_neighbours = self.getProperty("SmoothNeighbours").value
        ms.SmoothData(InputWorkspace=total_scatter_correction,
                      OutputWorkspace=total_scatter_correction,
                      NPoints=smooth_neighbours)
        ms.SmoothData(InputWorkspace=multi_scatter_correction,
                      OutputWorkspace=multi_scatter_correction,
                      NPoints=smooth_neighbours)

        return total_scatter_correction, multi_scatter_correction

    def _calculate_hydrogen_intensity(self, hydrogen, constraints):
        material_builder = MaterialBuilder()
        hydrogen_cross_section = hydrogen.totalScatterXSection()
        hydrogen_intensity = 0
        default_weight = 1.0 / len(constraints)

        for symbol in constraints.keys():
            constraint = constraints[symbol].value
            material = material_builder.setFormula(symbol).build()
            cross_section = material.totalScatterXSection()
            cross_section_ratio = hydrogen_cross_section / cross_section
            weight = constraint.get('weight', default_weight).value
            factor = constraint.get('factor', 1).value
            hydrogen_intensity += cross_section_ratio * factor * weight * constraint['intensity'].value

        return hydrogen_intensity


# -----------------------------------------------------------------------------------------
AlgorithmFactory.subscribe(VesuvioCorrections)