Unverified Commit 8afa1cad authored by Jose Borreguero's avatar Jose Borreguero Committed by GitHub
Browse files

Merge pull request #30849 from mantidproject/30737_topaz_satellite_peaks_v3

TOPAZ: Problem with peak intensity assignment for satellite peaks
parents c35a74fe 5686ea06
...@@ -52,6 +52,7 @@ set(SRC_FILES ...@@ -52,6 +52,7 @@ set(SRC_FILES
src/ImportMDHistoWorkspace.cpp src/ImportMDHistoWorkspace.cpp
src/ImportMDHistoWorkspaceBase.cpp src/ImportMDHistoWorkspaceBase.cpp
src/Integrate3DEvents.cpp src/Integrate3DEvents.cpp
src/IntegrateQLabEvents.cpp
src/IntegrateEllipsoids.cpp src/IntegrateEllipsoids.cpp
src/IntegrateEllipsoidsTwoStep.cpp src/IntegrateEllipsoidsTwoStep.cpp
src/IntegrateFlux.cpp src/IntegrateFlux.cpp
...@@ -165,6 +166,7 @@ set( ...@@ -165,6 +166,7 @@ set(
inc/MantidMDAlgorithms/ImportMDHistoWorkspace.h inc/MantidMDAlgorithms/ImportMDHistoWorkspace.h
inc/MantidMDAlgorithms/ImportMDHistoWorkspaceBase.h inc/MantidMDAlgorithms/ImportMDHistoWorkspaceBase.h
inc/MantidMDAlgorithms/Integrate3DEvents.h inc/MantidMDAlgorithms/Integrate3DEvents.h
inc/MantidMDAlgorithms/IntegrateQLabEvents.h
inc/MantidMDAlgorithms/IntegrateEllipsoids.h inc/MantidMDAlgorithms/IntegrateEllipsoids.h
inc/MantidMDAlgorithms/IntegrateEllipsoidsTwoStep.h inc/MantidMDAlgorithms/IntegrateEllipsoidsTwoStep.h
inc/MantidMDAlgorithms/IntegrateFlux.h inc/MantidMDAlgorithms/IntegrateFlux.h
...@@ -275,6 +277,7 @@ set(TEST_FILES ...@@ -275,6 +277,7 @@ set(TEST_FILES
ImportMDEventWorkspaceTest.h ImportMDEventWorkspaceTest.h
ImportMDHistoWorkspaceTest.h ImportMDHistoWorkspaceTest.h
Integrate3DEventsTest.h Integrate3DEventsTest.h
IntegrateQLabEventsTest.h
IntegrateEllipsoidsTest.h IntegrateEllipsoidsTest.h
IntegrateEllipsoidsTwoStepTest.h IntegrateEllipsoidsTwoStepTest.h
IntegrateEllipsoidsWithSatellitesTest.h IntegrateEllipsoidsWithSatellitesTest.h
......
...@@ -11,7 +11,7 @@ ...@@ -11,7 +11,7 @@
#include "MantidDataObjects/EventWorkspace.h" #include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/PeaksWorkspace.h" #include "MantidDataObjects/PeaksWorkspace.h"
#include "MantidDataObjects/Workspace2D.h" #include "MantidDataObjects/Workspace2D.h"
#include "MantidMDAlgorithms/Integrate3DEvents.h" #include "MantidMDAlgorithms/IntegrateQLabEvents.h"
#include "MantidMDAlgorithms/MDTransfInterface.h" #include "MantidMDAlgorithms/MDTransfInterface.h"
#include "MantidMDAlgorithms/MDWSDescription.h" #include "MantidMDAlgorithms/MDWSDescription.h"
#include "MantidMDAlgorithms/UnitsConversionHelper.h" #include "MantidMDAlgorithms/UnitsConversionHelper.h"
...@@ -24,30 +24,43 @@ namespace MDAlgorithms { ...@@ -24,30 +24,43 @@ namespace MDAlgorithms {
class DLLExport IntegrateEllipsoids : public API::Algorithm { class DLLExport IntegrateEllipsoids : public API::Algorithm {
public: public:
const std::string name() const override; const std::string name() const override { return "IntegrateEllipsoids"; }
/// Summary of algorithms purpose
const std::string summary() const override { const std::string summary() const override {
return "Integrate Single Crystal Diffraction Bragg peaks using 3D " return "Integrate Single Crystal Diffraction Bragg peaks using 3D "
"ellipsoids."; "ellipsoids.";
} }
int version() const override; int version() const override { return 1; }
const std::vector<std::string> seeAlso() const override { const std::vector<std::string> seeAlso() const override {
return {"IntegrateEllipsoidsTwoStep"}; return {"IntegrateEllipsoidsTwoStep"};
} }
const std::string category() const override; const std::string category() const override { return "Crystal\\Integration"; }
private: private:
using PrincipleAxes = std::array<std::vector<double>, 3>; /// Initialize the algorithm's properties
void init() override; void init() override;
/// Execute the algorithm
void exec() override; void exec() override;
void qListFromEventWS(Integrate3DEvents &integrator, API::Progress &prog,
DataObjects::EventWorkspace_sptr &wksp, /**
Kernel::DblMatrix const &UBinv, bool hkl_integ); * @brief create a list of SlimEvent objects from an events workspace
void qListFromHistoWS(Integrate3DEvents &integrator, API::Progress &prog, * @param integrator : integrator object on the list is accumulated
DataObjects::Workspace2D_sptr &wksp, * @param prog : progress object
Kernel::DblMatrix const &UBinv, bool hkl_integ); * @param wksp : input EventWorkspace
*/
void qListFromEventWS(IntegrateQLabEvents &integrator, API::Progress &prog,
DataObjects::EventWorkspace_sptr &wksp);
/**
* @brief create a list of SlimEvent objects from a histogram workspace
* @param integrator : integrator object on which the list is accumulated
* @param prog : progress object
* @param wksp : input Workspace2D
*/
void qListFromHistoWS(IntegrateQLabEvents &integrator, API::Progress &prog,
DataObjects::Workspace2D_sptr &wksp);
/// Calculate if this Q is on a detector /// Calculate if this Q is on a detector
void calculateE1(const Geometry::DetectorInfo &detectorInfo); void calculateE1(const Geometry::DetectorInfo &detectorInfo);
...@@ -60,6 +73,10 @@ private: ...@@ -60,6 +73,10 @@ private:
MDWSDescription m_targWSDescr; MDWSDescription m_targWSDescr;
/**
* @brief Initialize the output information for the MD conversion framework.
* @param wksp : The workspace to get information from.
*/
void initTargetWSDescr(API::MatrixWorkspace_sptr &wksp); void initTargetWSDescr(API::MatrixWorkspace_sptr &wksp);
}; };
......
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 2012 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 +
#pragma once
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/Peak.h"
#include "MantidDataObjects/PeakShapeEllipsoid.h"
#include "MantidKernel/Matrix.h"
#include "MantidKernel/V3D.h"
#include <memory>
#include <tuple>
#include <unordered_map>
#include <vector>
namespace Mantid {
namespace Geometry {
class PeakShape;
}
namespace MDAlgorithms {
using Mantid::DataObjects::PeakShapeEllipsoid_const_sptr;
using Mantid::Geometry::PeakShape_const_sptr;
using Mantid::Kernel::V3D;
/// Partition QLab space into a cubic lattice
struct CellCoords {
int64_t a;
int64_t b;
int64_t c;
CellCoords(const V3D &q, const double cellSize)
: a(static_cast<int64_t>(q[0] / cellSize)),
b(static_cast<int64_t>(q[1] / cellSize)),
c(static_cast<int64_t>(q[2] / cellSize)) {}
/// Check if all cell coords are zero
bool isOrigin() { return !(a || b || c); }
/// cast coordinates to scalar, to be used as key in an unordered map
int64_t getHash() { return 1000000000000 * a + 100000000 * b + 10000 * c; }
/// Hashes for the 26 first neighbor coordinates plus the coordinates
/// themselves
std::vector<int64_t> nearbyCellHashes() {
std::vector<int64_t> neighbors;
for (int64_t ia = a - 1; ia <= a + 1; ia++)
for (int64_t ib = b - 1; ib <= b + 1; ib++)
for (int64_t ic = c - 1; ic <= c + 1; ic++) {
int64_t key = 1000000000000 * ia + 100000000 * ib + 10000 * ic;
neighbors.emplace_back(key);
}
return neighbors;
}
};
// [(weight, error), Q-vector], trimmed-down info for an event
using SlimEvent = std::pair<std::pair<double, double>, V3D>;
using SlimEvents = std::vector<SlimEvent>;
// A cell in partitioned QLab space containing one peak
struct OccupiedCell {
size_t peakIndex; // index of the peak within this cell
V3D peakQ; // QLab vector of the peak within this cell
SlimEvents events; // events potentially closer than m_radius to the peak
};
/**
@class IntegrateQLabEvents
This is a low-level class to construct a map with lists of events near
each peak Q-vector in the lab frame. The Q-vector of each event is shifted
by the Q-vector of the associated peak.
A method is also provided to find the principal axes of such a list
of events and to find the net integrated counts using ellipsoids
with axis lengths determined from the standard deviations in the
directions of the principal axes.
@author Dennis Mikkelson
@date 2012-12-19
*/
class DLLExport IntegrateQLabEvents {
public:
/**
* @brief Store events within a certain radius of the specified peak centers,
* and sum these events to estimate pixel intensities.
* @param peak_q_list : List of Q-vectors for peak centers.
* @param radius : The maximum distance from a peak's Q-vector, for an
* event to be stored in the list associated with that peak.
* @param useOnePercentBackgroundCorrection : flag if one percent background
* correction should be used. */
IntegrateQLabEvents(const SlimEvents &peak_q_list, double radius,
const bool useOnePercentBackgroundCorrection = true);
/// Determine if an input Q-vector lies in the cell associated to the origin
static bool isOrigin(const V3D &q, const double &cellSize);
/**
* @brief distribute the events among the cells of the partitioned QLab
* space.
* @details Given QLab partitioned into a cubic lattice with unit cell of
* certain size, assign each event to one particular cell depending on its
* QLab vector.
* @param event_qs : List of SlimEvent objects to be distributed */
void addEvents(SlimEvents const &event_qs);
/**
* @brief Integrate the events around the specified peak QLab vector.
* @details The principal axes of the events near this Q-vector
* and the standard deviations in the directions of these principal
* axes determine ellipsoidal regions for integrating the peak and
* estimating the background. Alternatively, if peak and background
* radii are specified, then those will be used for half the major
* axis length of the ellipsoids, and the other axes of the ellipsoids
* will be set proportionally, based on the standard deviations.
* @param E1Vec : Vector of values for calculating edge of detectors
* @param peak_q : The QLab-vector for the peak center.
* @param specify_size : If true the integration will be done using the
* ellipsoids with major axes determined by the peak, back_inner
* and back_outer radii parameters. If false, the integration will be
* done using a peak region with major axis chosen so that it
* covers +- three standard deviations of the data in each direction. In
* this case, the background ellipsoidal shell is chosen to have the
* same VOLUME as the peak ellipsoid, and to use the peak ellipsoid
* for the inner radius.
* @param peak_radius : Size of half the major axis of the ellipsoidal
* peak region.
* @param back_inner_radius : Size of half the major axis of the INNER
* ellipsoidal boundary of the background region
* @param back_outer_radius : Size of half the major axis of the OUTER
* ellipsoidal boundary of the background region
* @param axes_radii : The radii used for integration in the directions
* of the three principal axes.
* @param inti : (output) collects the net integrated intensity
* @param sigi : (output) collects an estimate of the standard deviation
* of the net integrated intensity */
PeakShape_const_sptr ellipseIntegrateEvents(
const std::vector<V3D> &E1Vec, V3D const &peak_q, bool specify_size,
double peak_radius, double back_inner_radius, double back_outer_radius,
std::vector<double> &axes_radii, double &inti, double &sigi);
/**
* @brief Assign events to each of the cells occupied by events.
* @details Iterate over each QLab cell containing a peak and accumulate the
* list of events for the cell and for the first-neighbor cells into a
* single list of events. The QLab vectors for this events are shifted
* by the QLab vector of the peak. */
void populateCellsWithPeaks();
private:
/**
* @brief Number of events in an ellipsoid.
* @details The ellipsoid is centered at 0,0,0 with the three specified
* axes and the three specified sizes in the direction of those axes.
* NOTE: The three axes must be mutually orthogonal unit vectors.
* @param events : List of SlimEvents centered at 0,0,0
* @param directions : List of 3 orthonormal directions for the axes of
* the ellipsoid.
* @param sizes : List of three values a,b,c giving half the length
* of the three axes of the ellisoid.
* @return number of events and estimated error */
static std::pair<double, double>
numInEllipsoid(SlimEvents const &events, std::vector<V3D> const &directions,
std::vector<double> const &sizes);
/**
* @brief Number of events in an ellipsoid with background correction.
* @details The ellipsoid is centered at 0,0,0 with the three specified
* axes and the three specified sizes in the direction of those axes.
* NOTE: The three axes must be mutually orthogonal unit vectors.
* @param events : List of 3D events centered at 0,0,0
* @param directions : List of 3 orthonormal directions for the axes of
* the ellipsoid.
* @param sizes : List of three values a,b,c giving half the length
* of the three axes of the ellisoid.
* @param sizesIn : List of three values a,b,c giving half the length
* of the three inner axes of the ellisoid.
* @param useOnePercentBackgroundCorrection : flag if one percent background
* correction should be used.
* @return number of events and estimated error */
static std::pair<double, double> numInEllipsoidBkg(
SlimEvents const &events, std::vector<V3D> const &directions,
std::vector<double> const &sizes, std::vector<double> const &sizesIn,
const bool useOnePercentBackgroundCorrection);
/**
* @brief 3x3 covariance matrix of a list of SlimEvent objects
* @details the purpose of the covariance matrix is to find the principal axes
* of the SlimeEvents, associated with a particular peak. Their QLab vectors
* are already shifted by the QLab vector of the peak. Only events within
* the specified distance from the peak (here at Q=[0,0,0]) will be used.
* The covariance matrix can be easily constructed. X, Y, Z of each peak
* position are the variables we wish to determine the covariance. The mean
* position in each dimension has already been calculated on subtracted,
* since this corresponds to the QLab vector peak. The expected values
* of each correlation test X,X X,Y X,Z e.t.c form the elements of this
* 3x3 matrix, but since the probabilities are equal, we can remove them
* from the sums of the expected values, and simply divide by the number
* of events for each matrix element. Note that the diagonal elements
* form the variance X,X, Y,Y, Z,Z
* @param events : SlimEvents associated to one peak
* @param matrix : (output) 3x3 covariance matrix
* @param radius : Only events within this distance radius of the
* peak (here at Q=[0,0,0]) are used for calculating the covariance matrix.*/
static void makeCovarianceMatrix(SlimEvents const &events,
Kernel::DblMatrix &matrix, double radius);
/**
* @brief Eigen vectors of a 3x3 real symmetric matrix using the GSL.
* @param cov_matrix : 3x3 real symmetric matrix.
* @param eigen_vectors : (output) returned eigen vectors
* @param eigen_values : (output) three eigenvalues
*/
static void getEigenVectors(Kernel::DblMatrix const &cov_matrix,
std::vector<V3D> &eigen_vectors,
std::vector<double> &eigen_values);
/**
* @brief assign an event to one cell of the partitioned QLab space.
* @param event : SlimEvent to be assigned */
void addEvent(const SlimEvent event);
/**
* @brief Integrate a list of events associated to one peak.
* @details The QLab vector of the events are shifted by the QLab vector
* of the peak. Spatial distribution of the events in QLab space is
* described with principal axes of the ellipsoid, as well as the
* standard deviations in the the directions of the principal axes.
* @param E1Vec : Vector of values for calculating edge of detectors
* @param peak_q : The Q-vector for the peak center.
* @param ev_list : List of events centered around the peak (here with
* Q=[0,0,0]).
* @param directions : The three principal axes of the list of events
* @param sigmas : The standard deviations of the events in the
* directions of the three principal axes.
* @param specify_size : If true the integration will be done using the
* ellipsoids with major axes determined by the peak, back_inner and
* back_outer radii parameters. If false, the integration will be done
* using a peak region with major axis chosen so that it covers +- three
* standard deviations of the data in each direction. In this case, the
* background ellipsoidal shell is chosen to have the same VOLUME as the
* peak ellipsoid, and to use the peak ellipsoid for the inner radius.
* @param peak_radius : Size of half the major axis of the ellipsoid
* @param back_inner_radius : Size of half the major axis of the INNER
* ellipsoidal boundary of the background region
* @param back_outer_radius : Size of half the major axis of the OUTER
* ellipsoidal boundary of the background region
* @param axes_radii : The radii used for integration in the directions
* of the three principal axes.
* @param inti : (output) net integrated intensity
* @param sigi : (output) estimate of the standard deviation the intensity */
PeakShapeEllipsoid_const_sptr ellipseIntegrateEvents(
const std::vector<V3D> &E1Vec, V3D const &peak_q,
SlimEvents const &ev_list, std::vector<V3D> const &directions,
std::vector<double> const &sigmas, bool specify_size, double peak_radius,
double back_inner_radius, double back_outer_radius,
std::vector<double> &axes_radii, double &inti, double &sigi);
/**
* @brief Calculate if this Q is on a detector
* @details The distance from C to OE is given by dv=C-E*(C.scalar_prod(E))
* If dv.norm<integration_radius, one of the detector trajectories on the
* edge is too close to the peak. This method is applied to all masked
* pixels. If there are masked pixels trajectories inside an integration
* volume, the peak must be rejected.
* @param E1Vec : Vector of values for calculating edge of detectors
* @param QLabFrame: The Peak center.
* @param r: Peak radius.
*/
double detectorQ(const std::vector<V3D> &E1Vec, const V3D QLabFrame,
const std::vector<double> &r);
// Private data members
double m_radius; // size of sphere to use for events around a peak
/// if one percent culling of the background should be performed.
const bool m_useOnePercentBackgroundCorrection;
/// size of the square cell unit, holding at most one single peak
double m_cellSize;
/// list the occupied cells in an unordered map for fast searching
std::unordered_map<size_t, OccupiedCell> m_cellsWithPeaks;
/// list of cells occupied with events
std::unordered_map<size_t, SlimEvents> m_cellsWithEvents;
};
} // namespace MDAlgorithms
} // namespace Mantid
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 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 +
#include "MantidMDAlgorithms/IntegrateQLabEvents.h"
#include "MantidDataObjects/NoShape.h"
#include "MantidDataObjects/PeakShapeEllipsoid.h"
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidKernel/MultiThreaded.h"
#include <boost/math/special_functions/round.hpp>
#include <cmath>
#include <fstream>
#include <memory>
#include <numeric>
#include <tuple>
extern "C" {
#include <cstdio>
#include <utility>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
}
using namespace Mantid::DataObjects;
namespace Mantid {
using Mantid::API::Progress;
namespace MDAlgorithms {
using namespace std;
using Mantid::Kernel::DblMatrix;
using Mantid::Kernel::V3D;
IntegrateQLabEvents::IntegrateQLabEvents(
const SlimEvents &peak_q_list, double radius,
const bool useOnePercentBackgroundCorrection)
: m_radius(radius),
m_useOnePercentBackgroundCorrection(useOnePercentBackgroundCorrection) {
m_cellSize = m_radius;
for (size_t peakIndex = 0; peakIndex != peak_q_list.size(); ++peakIndex) {
const V3D q = peak_q_list[peakIndex].second;
CellCoords abc(q, m_cellSize);
// abc = [0, 0, 0] is no scattering
if (!abc.isOrigin()) {
SlimEvents events; // empty list
OccupiedCell c = {peakIndex, q, events};
std::pair<size_t, OccupiedCell> newCell(abc.getHash(), c);
m_cellsWithPeaks.insert(newCell);
}
}
}
bool IntegrateQLabEvents::isOrigin(const V3D &q, const double &cellSize) {
int64_t a(static_cast<int64_t>(q[0] / cellSize));
int64_t b(static_cast<int64_t>(q[1] / cellSize));
int64_t c(static_cast<int64_t>(q[2] / cellSize));
return !(a || b || c);
}
void IntegrateQLabEvents::addEvents(SlimEvents const &event_qs) {
for (const auto &event_q : event_qs)
addEvent(event_q);
}
Mantid::Geometry::PeakShape_const_sptr
IntegrateQLabEvents::ellipseIntegrateEvents(
const std::vector<V3D> &E1Vec, V3D const &peak_q, bool specify_size,
double peak_radius, double back_inner_radius, double back_outer_radius,
std::vector<double> &axes_radii, double &inti, double &sigi) {
inti = 0.0; // default values, in case something
sigi = 0.0; // is wrong with the peak.
int64_t hash = CellCoords(peak_q, m_cellSize).getHash();
auto cell_it = m_cellsWithPeaks.find(hash);
if (cell_it == m_cellsWithPeaks.end())
return std::make_shared<NoShape>(); // peak_q is [0, 0, 0]
OccupiedCell cell = cell_it->second;
SlimEvents &some_events = cell.events;
if (some_events.size() < 3)
return std::make_shared<NoShape>();
DblMatrix cov_matrix(3, 3);
makeCovarianceMatrix(some_events, cov_matrix, m_radius);
std::vector<V3D> eigen_vectors;
std::vector<double> eigen_values;
getEigenVectors(cov_matrix, eigen_vectors, eigen_values);
std::vector<double> sigmas(3);
for (int i = 0; i < 3; i++)
sigmas[i] = sqrt(eigen_values[i]);
bool invalid_peak =
std::any_of(sigmas.cbegin(), sigmas.cend(), [](const double sigma) {
return std::isnan(sigma) || sigma <= 0;
});
// if data collapses to a line or plane, the ellipsoid volume is zero
if (invalid_peak)
return std::make_shared<NoShape>();
return ellipseIntegrateEvents(std::move(E1Vec), peak_q, some_events,
eigen_vectors, sigmas, specify_size,
peak_radius, back_inner_radius,
back_outer_radius, axes_radii, inti, sigi);
}
std::pair<double, double>
IntegrateQLabEvents::numInEllipsoid(SlimEvents const &events,
std::vector<V3D> const &directions,
std::vector<double> const &sizes) {
std::pair<double, double> count(0, 0);
for (const auto &event : events) {
double sum = 0;
for (size_t k = 0; k < 3; k++) {
double comp = event.second.scalar_prod(directions[k]) / sizes[k];
sum += comp * comp;
}
if (sum <= 1) {
count.first += event.first.first; // count
count.second += event.first.second; // error squared (add in quadrature)
}
}
return count;
}
std::pair<double, double> IntegrateQLabEvents::numInEllipsoidBkg(
SlimEvents const &events, std::vector<V3D> const &directions,
std::vector<double> const &sizes, std::vector<double> const &sizesIn,
const bool useOnePercentBackgroundCorrection) {
std::pair<double, double> count(0, 0);
std::vector<std::pair<double, double>> eventVec;
for (const auto &event : events) {
double sum = 0;
double sumIn = 0;
for (size_t k = 0; k < 3; k++) {
double comp = event.second.scalar_prod(directions[k]) / sizes[k];
sum += comp * comp;
comp = event.second.scalar_prod(directions[k]) / sizesIn[k];
sumIn += comp * comp;
}
if (sum <= 1 && sumIn >= 1)
eventVec.emplace_back(event.first);
}
auto endIndex = eventVec.size();
if (useOnePercentBackgroundCorrection) {
// Remove top 1% of background
std::sort(
eventVec.begin(), eventVec.end(),
[](const std::pair<double, double> &a,
const std::pair<double, double> &b) { return a.first < b.first; });
endIndex = static_cast<size_t>(0.99 * static_cast<double>(endIndex));
}
for (size_t k = 0; k < endIndex; ++k) {
count.first += eventVec[k].first;
count.second += eventVec[k].second;
}