Skip to content
Snippets Groups Projects
Commit 0e11df64 authored by Savici, Andrei T.'s avatar Savici, Andrei T.
Browse files

MDNormDirectSC can now handle merged MDEventWorkspaces

parent bddad0bb
No related branches found
No related tags found
No related merge requests found
...@@ -55,13 +55,13 @@ private: ...@@ -55,13 +55,13 @@ private:
DataObjects::MDHistoWorkspace_sptr binInputWS(); DataObjects::MDHistoWorkspace_sptr binInputWS();
void createNormalizationWS(const DataObjects::MDHistoWorkspace &dataWS); void createNormalizationWS(const DataObjects::MDHistoWorkspace &dataWS);
std::vector<coord_t> std::vector<coord_t>
getValuesFromOtherDimensions(bool &skipNormalization) const; getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex=0) const;
Kernel::Matrix<coord_t> Kernel::Matrix<coord_t>
findIntergratedDimensions(const std::vector<coord_t> &otherDimValues, findIntergratedDimensions(const std::vector<coord_t> &otherDimValues,
bool &skipNormalization); bool &skipNormalization);
void cacheDimensionXValues(); void cacheDimensionXValues();
void calculateNormalization(const std::vector<coord_t> &otherValues, void calculateNormalization(const std::vector<coord_t> &otherValues,
const Kernel::Matrix<coord_t> &affineTrans); const Kernel::Matrix<coord_t> &affineTrans, uint16_t expInfoIndex);
void calculateIntersections(std::vector<std::array<double, 4>> &intersections, void calculateIntersections(std::vector<std::array<double, 4>> &intersections,
const double theta, const double phi); const double theta, const double phi);
...@@ -88,7 +88,10 @@ private: ...@@ -88,7 +88,10 @@ private:
Kernel::V3D m_beamDir; Kernel::V3D m_beamDir;
/// ki-kf for Inelastic convention; kf-ki for Crystallography convention /// ki-kf for Inelastic convention; kf-ki for Crystallography convention
std::string convention; std::string convention;
/// internal flag to accumulate to an existing workspace
bool m_accumulate{false}; bool m_accumulate{false};
/// number of experiment infos
uint16_t m_numExptInfos;
}; };
} // namespace MDAlgorithms } // namespace MDAlgorithms
......
...@@ -142,22 +142,31 @@ void MDNormDirectSC::exec() { ...@@ -142,22 +142,31 @@ void MDNormDirectSC::exec() {
m_normWS->setDisplayNormalization(Mantid::API::NoNormalization); m_normWS->setDisplayNormalization(Mantid::API::NoNormalization);
setProperty("OutputNormalizationWorkspace", m_normWS); setProperty("OutputNormalizationWorkspace", m_normWS);
// Check for other dimensions if we could measure anything in the original m_numExptInfos = outputWS->getNumExperimentInfo();
// data // loop over all experiment infos
bool skipNormalization = false; for (uint16_t expInfoIndex=0; expInfoIndex < m_numExptInfos; expInfoIndex++)
const std::vector<coord_t> otherValues = {
getValuesFromOtherDimensions(skipNormalization); if (expInfoIndex > 0) {
const auto affineTrans = m_accumulate = true;
findIntergratedDimensions(otherValues, skipNormalization); }
cacheDimensionXValues(); // Check for other dimensions if we could measure anything in the original
// data
if (!skipNormalization) { bool skipNormalization = false;
calculateNormalization(otherValues, affineTrans); const std::vector<coord_t> otherValues =
} else { getValuesFromOtherDimensions(skipNormalization,expInfoIndex);
g_log.warning("Binning limits are outside the limits of the MDWorkspace. " const auto affineTrans =
"Not applying normalization."); findIntergratedDimensions(otherValues, skipNormalization);
cacheDimensionXValues();
if (!skipNormalization) {
calculateNormalization(otherValues, affineTrans,expInfoIndex);
} else {
g_log.warning("Binning limits are outside the limits of the MDWorkspace. "
"Not applying normalization.");
}
} }
// Set the display normalization based on the input workspace // Set the display normalization based on the input workspace
outputWS->setDisplayNormalization(m_inputWS->displayNormalizationHisto()); outputWS->setDisplayNormalization(m_inputWS->displayNormalizationHisto());
} }
...@@ -299,8 +308,8 @@ void MDNormDirectSC::createNormalizationWS(const MDHistoWorkspace &dataWS) { ...@@ -299,8 +308,8 @@ void MDNormDirectSC::createNormalizationWS(const MDHistoWorkspace &dataWS) {
* MD position calculation * MD position calculation
*/ */
std::vector<coord_t> std::vector<coord_t>
MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization) const { MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization, uint16_t expInfoIndex) const {
const auto &runZero = m_inputWS->getExperimentInfo(0)->run(); const auto &currentRun = m_inputWS->getExperimentInfo(expInfoIndex)->run();
std::vector<coord_t> otherDimValues; std::vector<coord_t> otherDimValues;
for (size_t i = 4; i < m_inputWS->getNumDims(); i++) { for (size_t i = 4; i < m_inputWS->getNumDims(); i++) {
...@@ -308,7 +317,7 @@ MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization) const { ...@@ -308,7 +317,7 @@ MDNormDirectSC::getValuesFromOtherDimensions(bool &skipNormalization) const {
float dimMin = static_cast<float>(dimension->getMinimum()); float dimMin = static_cast<float>(dimension->getMinimum());
float dimMax = static_cast<float>(dimension->getMaximum()); float dimMax = static_cast<float>(dimension->getMaximum());
auto *dimProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>( auto *dimProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
runZero.getProperty(dimension->getName())); currentRun.getProperty(dimension->getName()));
if (dimProp) { if (dimProp) {
coord_t value = static_cast<coord_t>(dimProp->firstValue()); coord_t value = static_cast<coord_t>(dimProp->firstValue());
otherDimValues.push_back(value); otherDimValues.push_back(value);
...@@ -396,7 +405,8 @@ Kernel::Matrix<coord_t> MDNormDirectSC::findIntergratedDimensions( ...@@ -396,7 +405,8 @@ Kernel::Matrix<coord_t> MDNormDirectSC::findIntergratedDimensions(
} }
/** /**
* Stores the X values from each H,K,L dimension as member variables * Stores the X values from each H,K,L,E dimension as member variables
* Energy dimension is transformed to final wavevector.
*/ */
void MDNormDirectSC::cacheDimensionXValues() { void MDNormDirectSC::cacheDimensionXValues() {
constexpr double energyToK = 8.0 * M_PI * M_PI * constexpr double energyToK = 8.0 * M_PI * M_PI *
...@@ -442,17 +452,17 @@ void MDNormDirectSC::cacheDimensionXValues() { ...@@ -442,17 +452,17 @@ void MDNormDirectSC::cacheDimensionXValues() {
* @param otherValues * @param otherValues
* @param affineTrans * @param affineTrans
*/ */
void MDNormDirectSC::calculateNormalization( void MDNormDirectSC::calculateNormalization(const std::vector<coord_t> &otherValues,
const std::vector<coord_t> &otherValues, const Kernel::Matrix<coord_t> &affineTrans,
const Kernel::Matrix<coord_t> &affineTrans) { uint16_t expInfoIndex) {
constexpr double energyToK = 8.0 * M_PI * M_PI * constexpr double energyToK = 8.0 * M_PI * M_PI *
PhysicalConstants::NeutronMass * PhysicalConstants::NeutronMass *
PhysicalConstants::meV * 1e-20 / PhysicalConstants::meV * 1e-20 /
(PhysicalConstants::h * PhysicalConstants::h); (PhysicalConstants::h * PhysicalConstants::h);
const auto &exptInfoZero = *(m_inputWS->getExperimentInfo(0)); const auto &currentExptInfo = *(m_inputWS->getExperimentInfo(expInfoIndex));
using VectorDoubleProperty = Kernel::PropertyWithValue<std::vector<double>>; using VectorDoubleProperty = Kernel::PropertyWithValue<std::vector<double>>;
auto *rubwLog = auto *rubwLog =
dynamic_cast<VectorDoubleProperty *>(exptInfoZero.getLog("RUBW_MATRIX")); dynamic_cast<VectorDoubleProperty *>(currentExptInfo.getLog("RUBW_MATRIX"));
if (!rubwLog) { if (!rubwLog) {
throw std::runtime_error( throw std::runtime_error(
"Wokspace does not contain a log entry for the RUBW matrix." "Wokspace does not contain a log entry for the RUBW matrix."
...@@ -460,12 +470,12 @@ void MDNormDirectSC::calculateNormalization( ...@@ -460,12 +470,12 @@ void MDNormDirectSC::calculateNormalization(
} else { } else {
Kernel::DblMatrix rubwValue( Kernel::DblMatrix rubwValue(
(*rubwLog)()); // includes the 2*pi factor but not goniometer for now :) (*rubwLog)()); // includes the 2*pi factor but not goniometer for now :)
m_rubw = exptInfoZero.run().getGoniometerMatrix() * rubwValue; m_rubw = currentExptInfo.run().getGoniometerMatrix() * rubwValue;
m_rubw.Invert(); m_rubw.Invert();
} }
const double protonCharge = exptInfoZero.run().getProtonCharge(); const double protonCharge = currentExptInfo.run().getProtonCharge();
const auto &spectrumInfo = exptInfoZero.spectrumInfo(); const auto &spectrumInfo = currentExptInfo.spectrumInfo();
// Mapping // Mapping
const int64_t ndets = static_cast<int64_t>(spectrumInfo.size()); const int64_t ndets = static_cast<int64_t>(spectrumInfo.size());
...@@ -482,7 +492,8 @@ void MDNormDirectSC::calculateNormalization( ...@@ -482,7 +492,8 @@ void MDNormDirectSC::calculateNormalization(
std::vector<std::atomic<signal_t>> signalArray(m_normWS->getNPoints()); std::vector<std::atomic<signal_t>> signalArray(m_normWS->getNPoints());
std::vector<std::array<double, 4>> intersections; std::vector<std::array<double, 4>> intersections;
std::vector<coord_t> pos, posNew; std::vector<coord_t> pos, posNew;
auto prog = make_unique<API::Progress>(this, 0.3, 1.0, ndets); double progStep = 0.7/m_numExptInfos;
auto prog = make_unique<API::Progress>(this, 0.3 + progStep * expInfoIndex, 0.3 + progStep * (expInfoIndex + 1.), ndets);
// cppcheck-suppress syntaxError // cppcheck-suppress syntaxError
PRAGMA_OMP(parallel for private(intersections, pos, posNew)) PRAGMA_OMP(parallel for private(intersections, pos, posNew))
for (int64_t i = 0; i < ndets; i++) { for (int64_t i = 0; i < ndets; i++) {
......
...@@ -14,17 +14,15 @@ The algorithm calculates a normalization MD workspace for single crystal direct ...@@ -14,17 +14,15 @@ The algorithm calculates a normalization MD workspace for single crystal direct
Trajectories of each detector in reciprocal space are calculated, and the flux is integrated between intersections with each Trajectories of each detector in reciprocal space are calculated, and the flux is integrated between intersections with each
MDBox. A brief introduction to the multi-dimensional data normalization can be found :ref:`here <MDNorm>`. MDBox. A brief introduction to the multi-dimensional data normalization can be found :ref:`here <MDNorm>`.
.. Note::
This is an experimental algorithm in Release 3.3. Please check the nightly Mantid build, and the Mantid webpage
for better offline help and usage examples.
.. Note:: .. Note::
If the MDEvent input workspace is generated from an event workspace, the algorithm gives the correct normalization If the MDEvent input workspace is generated from an event workspace, the algorithm gives the correct normalization
only if the event workspace is cropped and binned to the same energy transfer range. If the workspace is not cropped, only if the event workspace is cropped and binned to the same energy transfer range. If the workspace is not cropped,
one might have events in places where the normalization is calculated to be 0. one might have events in places where the normalization is calculated to be 0.
.. Note::
As of :ref:`Release 3.14.0 <v3.14.0>`, the algorithm can handle merged MD workspaces. Make sure all original MDEvent workspaces have the same dimensions
Usage Usage
----- -----
......
...@@ -25,11 +25,8 @@ Improvements ...@@ -25,11 +25,8 @@ Improvements
- Improved ``Save``-section of the TOFTOF reduction dialog. - Improved ``Save``-section of the TOFTOF reduction dialog.
- Behavior of the :ref:`LoadDNSLegacy <algm-LoadDNSLegacy>` for TOF data has been changed: the algorithm does not try to guess elastic channel any more, but asks for the user input. - Behavior of the :ref:`LoadDNSLegacy <algm-LoadDNSLegacy>` for TOF data has been changed: the algorithm does not try to guess elastic channel any more, but asks for the user input.
- :ref:`LoadDNSSCD <algm-LoadDNSSCD>` has been improved to be able to load TOF data.
- :ref:`MDNormDirectSC <algm-MDNormDirectSC>` now can handle merged MD workspaces.
:ref:`Release 3.14.0 <v3.14.0>` :ref:`Release 3.14.0 <v3.14.0>`
Improvements
############
- :ref:`LoadDNSSCD <algm-LoadDNSSCD>` has been improved to be able to load TOF data.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment