28#include <boost/math/special_functions/round.hpp>
53 auto ws_valid = std::make_shared<CompositeValidator>();
56 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
57 mustBePositive->setLower(0.0);
61 "An input MatrixWorkspace with time-of-flight units along "
62 "X-axis and defined instrument with defined sample");
65 "Workspace with peaks to be integrated");
68 "Only events at most this distance from a peak will be "
69 "considered when integrating");
71 declareProperty(
"SpecifySize",
false,
"If true, use the following for the major axis sizes, else use 3-sigma");
73 declareProperty(
"PeakSize", .18, mustBePositive,
"Half-length of major axis for peak ellipsoid");
76 "Half-length of major axis for inner ellipsoidal surface of "
80 "Half-length of major axis for outer ellipsoidal surface of "
83 declareProperty(
"IntegrateInHKL",
false,
"If true, integrate in HKL space not Q space.");
85 "Set to false to not integrate if peak radius is off edge of detector."
86 "Background will be scaled if background radius is off edge.");
89 "Default is false. If true, "
90 "BackgroundOuterRadius + AdaptiveQMultiplier * **|Q|** and "
91 "BackgroundInnerRadius + AdaptiveQMultiplier * **|Q|**");
94 "PeakRadius + AdaptiveQMultiplier * **|Q|** "
96 "different integration radius. Q includes the 2*pi factor.");
98 declareProperty(
"WeakPeakThreshold", 1.0, mustBePositive,
"Intensity threshold use to classify a peak as weak.");
101 "If this options is enabled, then the the top 1% of the "
102 "background will be removed"
103 "before the background subtraction.");
106 "The output PeaksWorkspace will be a copy of the input PeaksWorkspace "
107 "with the peaks' integrated intensities.");
115 Workspace2D_sptr histoWS = std::dynamic_pointer_cast<Workspace2D>(input_ws);
116 if (!eventWS && !histoWS) {
117 throw std::runtime_error(
"IntegrateEllipsoids needs either a "
118 "EventWorkspace or Workspace2D as input.");
121 const double weakPeakThreshold =
getProperty(
"WeakPeakThreshold");
124 if (!input_peak_ws) {
125 throw std::runtime_error(
"Could not read the Peaks Workspace");
129 throw std::runtime_error(
"Could not read the Input Workspace");
133 if (peak_ws != input_peak_ws) {
134 peak_ws = input_peak_ws->clone();
137 Progress prog(
this, 0.5, 1.0, input_ws->getNumberHistograms());
139 std::vector<Peak> &peaks = peak_ws->getPeaks();
140 size_t n_peaks = peak_ws->getNumberPeaks();
141 size_t indexed_count = 0;
142 std::vector<V3D> peak_q_list;
143 std::vector<V3D> hkl_vectors;
144 for (
size_t i = 0; i < n_peaks; i++)
146 V3D hkl(peaks[i].getH(), peaks[i].getK(), peaks[i].getL());
150 peak_q_list.emplace_back(peaks[i].getQLabFrame());
151 V3D miller_ind(
static_cast<double>(boost::math::iround<double>(hkl[0])),
152 static_cast<double>(boost::math::iround<double>(hkl[1])),
153 static_cast<double>(boost::math::iround<double>(hkl[2])));
154 hkl_vectors.emplace_back(miller_ind);
159 if (indexed_count < 3) {
160 throw std::runtime_error(
"At least three linearly independent indexed peaks are needed.");
168 UBinv *= (1.0 / (2.0 * M_PI));
170 std::vector<std::pair<std::pair<double, double>,
V3D>> qList;
171 for (
size_t i = 0; i < n_peaks; i++) {
172 qList.emplace_back(std::pair<double, double>(1.0, 1.0),
V3D(peaks[i].getQLabFrame()));
175 const bool integrateEdge =
getProperty(
"IntegrateIfOnEdge");
176 if (!integrateEdge) {
183 g_log.
error(
"Can't execute MaskBTP algorithm for this instrument to set "
184 "edge for IntegrateIfOnEdge option");
189 const bool integrateInHKL =
getProperty(
"IntegrateInHKL");
190 bool useOnePercentBackgroundCorrection =
getProperty(
"UseOnePercentBackgroundCorrection");
201 std::vector<std::pair<int, V3D>> weakPeaks, strongPeaks;
204 for (
int index = 0;
static_cast<size_t>(
index) < qList.size(); ++
index) {
205 const auto &item = qList[
index];
206 const auto center = item.second;
210 auto &peak = peak_ws->getPeak(
index);
211 peak.setIntensity(0);
212 peak.setSigmaIntensity(0);
214 const auto result = std::make_pair(
index, center);
215 if (sig2noise < weakPeakThreshold) {
216 g_log.
notice() <<
"Peak " << peak.getHKL() <<
" with Q = " << center <<
" is a weak peak with signal to noise "
217 << sig2noise <<
"\n";
218 weakPeaks.emplace_back(result);
220 g_log.
notice() <<
"Peak " << peak.getHKL() <<
" with Q = " << center <<
" is a strong peak with signal to noise "
221 << sig2noise <<
"\n";
222 strongPeaks.emplace_back(result);
226 std::vector<std::pair<std::shared_ptr<const Geometry::PeakShape>, std::tuple<double, double, double>>> shapeLibrary;
229 for (
const auto &item : strongPeaks) {
230 const auto index = item.first;
231 const auto &q = item.second;
236 shapeLibrary.emplace_back(result);
238 auto &peak = peak_ws->getPeak(
index);
239 peak.setIntensity(inti);
240 peak.setSigmaIntensity(sigi);
241 peak.setPeakShape(std::get<0>(result));
244 std::vector<Eigen::Vector3d> points;
245 std::transform(strongPeaks.begin(), strongPeaks.end(), std::back_inserter(points),
246 [&](
const std::pair<int, V3D> &item) {
247 const auto q = item.second;
248 return Eigen::Vector3d(q[0], q[1], q[2]);
252 throw std::runtime_error(
"Cannot integrate peaks when all peaks are below "
253 "the signal to noise ratio.");
258 for (
const auto &item : weakPeaks) {
260 const auto index = item.first;
261 const auto &q = item.second;
263 const auto result = kdTree.
findNearest(Eigen::Vector3d(q[0], q[1], q[2]));
264 const auto strongIndex =
static_cast<int>(std::get<1>(result[0]));
266 auto &peak = peak_ws->getPeak(
index);
267 auto &strongPeak = peak_ws->getPeak(strongIndex);
269 g_log.
notice() <<
"Integrating weak peak " << peak.getHKL() <<
" using strong peak " << strongPeak.getHKL() <<
"\n";
271 const auto libShape = shapeLibrary[
static_cast<int>(strongIndex)];
272 const auto shape = std::dynamic_pointer_cast<const PeakShapeEllipsoid>(libShape.first);
273 const auto frac = std::get<0>(libShape.second);
275 g_log.
notice() <<
"Weak peak will be adjusted by " << frac <<
"\n";
277 const auto weakShape = integrator.
integrateWeakPeak(params, shape, libShape.second, q, inti, sigi);
279 peak.setIntensity(inti);
280 peak.setSigmaIntensity(sigi);
281 peak.setPeakShape(weakShape);
286 peak_ws->mutableRun().addProperty(
"PeaksIntegrated", 1,
true);
299 const bool adaptiveQBackground =
getProperty(
"AdaptiveQBackground");
300 const double adaptiveQMultiplier =
getProperty(
"AdaptiveQMultiplier");
301 const double adaptiveQBackgroundMultiplier = (adaptiveQBackground) ? adaptiveQMultiplier : 0.0;
304 const double lenQpeak = peak_q.
norm();
316 const std::string
ELASTIC(
"Elastic");
318 const std::string
Q3D(
"Q3D");
319 const std::size_t
DIMS(3);
322 m_targWSDescr.
setMinMax(std::vector<double>(3, -2000.), std::vector<double>(3, 2000.));
329 childAlg->setProperty(
"InputWorkspace", wksp);
330 childAlg->executeAsChildAlg();
334 throw(std::runtime_error(
"Can not retrieve results of \"PreprocessDetectorsToMD\""));
338 auto numSpectra =
static_cast<int>(wksp->getNumberHistograms());
340 for (
int i = 0; i < numSpectra; ++i) {
345 unitConverter.
initialize(m_targWSDescr,
"Momentum");
351 std::vector<double> buffer(
DIMS);
353 EventList &events = wksp->getSpectrum(i);
359 if (events.
empty()) {
365 std::vector<Mantid::coord_t> locCoord(
DIMS, 0.);
366 unitConverter.updateConversion(i);
373 std::vector<std::pair<std::pair<double, double>,
V3D>> qList;
374 for (
const auto &raw_event : raw_events) {
375 double val = unitConverter.convertUnits(raw_event.tof());
377 for (
size_t dim = 0; dim <
DIMS; ++dim) {
378 buffer[dim] = locCoord[dim];
380 V3D qVec(buffer[0], buffer[1], buffer[2]);
383 qList.emplace_back(std::pair<double, double>(raw_event.m_weight, raw_event.m_errorSquared), qVec);
403 DblMatrix const &UBinv,
bool hkl_integ) {
406 const std::string
ELASTIC(
"Elastic");
408 const std::string
Q3D(
"Q3D");
409 const std::size_t
DIMS(3);
412 m_targWSDescr.
setMinMax(std::vector<double>(3, -2000.), std::vector<double>(3, 2000.));
419 childAlg->setProperty(
"InputWorkspace", wksp);
420 childAlg->executeAsChildAlg();
424 throw(std::runtime_error(
"Can not retrieve results of \"PreprocessDetectorsToMD\""));
428 auto numSpectra =
static_cast<int>(wksp->getNumberHistograms());
430 for (
int i = 0; i < numSpectra; ++i) {
435 unitConverter.
initialize(m_targWSDescr,
"Momentum");
442 const auto &xVals = wksp->points(i);
443 const auto &yVals = wksp->y(i);
444 const auto &eVals = wksp->e(i);
447 std::vector<Mantid::coord_t> locCoord(
DIMS, 0.);
448 unitConverter.updateConversion(i);
455 std::vector<std::pair<std::pair<double, double>,
V3D>> qList;
457 for (
size_t j = 0; j < yVals.size(); ++j) {
458 const double &yVal = yVals[j];
459 const double &esqVal = eVals[j] * eVals[j];
462 double val = unitConverter.convertUnits(xVals[j]);
464 V3D qVec(locCoord[0], locCoord[1], locCoord[2]);
468 if (std::isnan(qVec[0]) || std::isnan(qVec[1]) || std::isnan(qVec[2]))
472 qList.emplace_back(std::pair<double, double>(yVal, esqVal), qVec);
495 for (
size_t i = 0; i < detectorInfo.
size(); ++i) {
500 const auto &det = detectorInfo.
detector(i);
502 double ph1 = det.getPhi();
503 V3D E1 =
V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1),
505 E1 =
E1 * (1. /
E1.norm());
511 const std::string &property,
const std::string &values) {
514 alg->setProperty(property, values);
516 throw std::runtime_error(
"MaskDetectors Child Algorithm has not executed successfully");
#define DECLARE_ALGORITHM(classname)
std::map< DeltaEMode::Type, std::string > index
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_CRITICAL(name)
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
virtual std::shared_ptr< Algorithm > createChildAlgorithm(const std::string &name, const double startProgress=-1., const double endProgress=-1., const bool enableLogging=true, const int &version=-1)
Create a Child Algorithm.
A validator which checks that a workspace has a valid instrument.
Helper class for reporting progress from algorithms.
A property class for workspaces.
void compressEvents(double tolerance, EventList *destination)
Compress the event list by grouping events with the same TOF (within a given tolerance).
std::vector< WeightedEventNoTime > & getWeightedEventsNoTime()
Return the list of WeightedEvent contained.
void switchTo(Mantid::API::EventType newType) override
Switch the EventList to use the given EventType (TOF, WEIGHTED, or WEIGHTED_NOTIME)
bool empty() const
Much like stl containers, returns true if there is nothing in the event list.
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
bool isMasked(const size_t index) const
Returns true if the detector is masked.
const Geometry::IDetector & detector(const size_t index) const
Return a const reference to the detector with given index.
size_t size() const
Returns the size of the DetectorInfo, i.e., the number of detectors in the instrument.
bool isMonitor(const size_t index) const
Returns true if the detector is a monitor.
virtual double getTwoTheta(const Kernel::V3D &observer, const Kernel::V3D &axis) const =0
Gives the angle of this detector object with respect to an axis.
static double Optimize_UB(Kernel::DblMatrix &UB, const std::vector< Kernel::V3D > &hkl_vectors, const std::vector< Kernel::V3D > &q_vectors, std::vector< double > &sigabc)
Find the UB matrix that most nearly maps hkl to qxyz for 3 or more peaks.
static bool ValidIndex(const Kernel::V3D &hkl, double tolerance)
Check is hkl is within tolerance of integer (h,k,l) non-zero values.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
void error(const std::string &msg)
Logs at error level.
T Invert()
LU inversion routine.
NearestNeighbourResults findNearest(const VectorType &pos, const size_t k=1, const double error=0.0)
Find the k nearest neighbours to a given point.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
double norm() const noexcept
This is a low-level class to construct a map with lists of events near each peak Q-vector,...
std::pair< std::shared_ptr< const Mantid::Geometry::PeakShape >, std::tuple< double, double, double > > integrateStrongPeak(const IntegrationParameters ¶ms, const Kernel::V3D &peak_q, double &inti, double &sigi)
Find the net integrated intensity of a peak, using ellipsoidal volumes.
std::shared_ptr< const Geometry::PeakShape > integrateWeakPeak(const IntegrationParameters ¶ms, Mantid::DataObjects::PeakShapeEllipsoid_const_sptr shape, const std::tuple< double, double, double > &libPeak, const Mantid::Kernel::V3D &peak_q, double &inti, double &sigi)
double estimateSignalToNoiseRatio(const IntegrationParameters ¶ms, const Mantid::Kernel::V3D ¢er, bool forceSpherical=false, double sphericityTol=0.02)
void addEvents(std::vector< std::pair< std::pair< double, double >, Mantid::Kernel::V3D > > const &event_qs, bool hkl_integ)
Add event Q's to lists of events near peaks.
IntegrateEllipsoidsTwoStep provides a two pass peak integration algorithm.
void qListFromHistoWS(Integrate3DEvents &integrator, API::Progress &prog, DataObjects::Workspace2D_sptr &wksp, const Kernel::DblMatrix &UBinv, bool hkl_integ)
qListFromHistoWS creates qlist from input workspaces of type Workspace2D
void exec() override
Virtual method - must be overridden by concrete algorithm.
void init() override
Virtual method - must be overridden by concrete algorithm.
IntegrationParameters makeIntegrationParameters(const Kernel::V3D &peak_q) const
void qListFromEventWS(Integrate3DEvents &integrator, API::Progress &prog, DataObjects::EventWorkspace_sptr &wksp, const Kernel::DblMatrix &UBinv, bool hkl_integ)
int version() const override
Get the version of this algorithm.
const std::string category() const override
Get the category of this algorithm.
void runMaskDetectors(const Mantid::DataObjects::PeaksWorkspace_sptr &peakWS, const std::string &property, const std::string &values)
void calculateE1(const Geometry::DetectorInfo &detectorInfo)
Calculate if this Q is on a detector.
std::vector< Kernel::V3D > E1Vec
save all detector pixels
Class responsible for conversion of input workspace data into proper number of output dimensions for ...
bool calcMatrixCoord(const double &deltaEOrK0, std::vector< coord_t > &Coord, double &s, double &err) const override
Calculates 3D transformation of the variable coordinates and (if applicable) signal and error dependi...
void initialize(const MDWSDescription &ConvParams) override
function initalizes all variables necessary for converting workspace variables into MD variables in M...
bool calcYDepCoordinates(std::vector< coord_t > &Coord, size_t i) override
Method updates the value of preprocessed detector coordinates in Q-space, used by other functions.
helper class describes the properties of target MD workspace, which should be obtained as the result ...
void setMinMax(const std::vector< double > &minVal, const std::vector< double > &maxVal)
function sets up min-max values to the dimensions, described by the class
void buildFromMatrixWS(const API::MatrixWorkspace_sptr &pWS, const std::string &QMode, const std::string &dEMode, const std::vector< std::string > &dimPropertyNames=std::vector< std::string >())
method builds MD Event ws description from a matrix workspace and the transformations,...
DataObjects::TableWorkspace_const_sptr m_PreprDetTable
void setLorentsCorr(bool On=false)
do we need to perform Lorentz corrections
void initialize(const MDWSDescription &targetWSDescr, const std::string &unitsTo, bool forceViaTOF=false)
Initialize unit conversion helper This method is interface to internal initialize method,...
std::complex< double > MANTID_API_DLL E1(std::complex< double > z)
Integral for Gamma.
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
const std::string Q3D("Q3D")
Only convert to Q-vector.
const std::size_t DIMS(3)
Q-vector is always three dimensional.
const std::string ELASTIC("Elastic")
This only works for diffraction.
@ InOut
Both an input & output workspace.
@ Input
An input workspace.
@ Output
An output workspace.
double backgroundOuterRadius
double backgroundInnerRadius
std::vector< Kernel::V3D > E1Vectors