34template <
typename MDE,
size_t nd>
struct IsFullEvent : boost::false_type {};
37template <
size_t nd>
struct IsFullEvent<
MDEvent<nd>, nd> : boost::true_type {};
43template <
typename MDE,
size_t nd>
bool isFullMDEvent(
const boost::true_type & ) {
return true; }
49template <
typename MDE,
size_t nd>
bool isFullMDEvent(
const boost::false_type & ) {
return false; }
54template <
typename MDE,
size_t nd>
bool isFullMDEvent() {
return isFullMDEvent<MDE, nd>(IsFullEvent<MDE, nd>()); }
61template <
typename MDE,
size_t nd>
64 std::cerr <<
"Box has children\n";
65 addDetectors(peak, box, boost::true_type());
69 throw std::invalid_argument(
"FindPeaksMD::addDetectors - Unexpected Box "
70 "type, cannot retrieve events");
72 const auto &events = mdBox->getConstEvents();
73 auto itend = events.end();
74 for (
auto it = events.begin(); it != itend; ++it) {
81template <
typename MDE,
size_t nd>
83 const boost::false_type & ) {
84 throw std::runtime_error(
"FindPeaksMD - Workspace contains lean events, "
85 "cannot include detector information");
95 addDetectors(peak, box, IsFullEvent<MDE, nd>());
109 : peakWS(), peakRadiusSquared(), DensityThresholdFactor(0.0), m_maxPeaks(0), m_addDetectors(true),
110 m_densityScaleFactor(1e-6), prog(nullptr), inst(), m_runNumber(-1), dimType(), m_goniometer() {}
117 "An input MDEventWorkspace or MDHistoWorkspace with at least "
121 "Threshold distance for rejecting peaks that are found to be too close "
123 "This should be some multiple of the radius of a peak. Default: 0.1.");
126 "Maximum number of peaks to find. Default: 500.");
130 "Strategy for finding peaks in an MD workspace."
131 "1. VolumeNormalization: This is the default strategy. It will sort "
132 "all boxes in the workspace by deacresing order of signal density "
133 "(total weighted event sum divided by box volume).\n"
134 "2.NumberOfEventsNormalization: This option is only valid for "
135 "MDEventWorkspaces. It will use the total weighted event sum divided"
136 "by the number of events. This can improve peak finding for "
138 "raw data which has been converted to an EventWorkspace. The threshold "
140 "peak finding can be controlled by the SingalThresholdFactor property "
142 "be larger than 1. Note that this approach does not work for event-based "
146 "The overall signal density of the workspace will be "
147 "multiplied by this factor \n"
148 "to get a threshold signal density below which boxes are NOT "
149 "considered to be peaks. See the help.\n"
153 std::make_unique<EnabledWhenProperty>(
157 "The overal signal value (not density!) normalized by the "
158 "number of events is compared to the specified signal "
159 "threshold. Boxes which are below this threshold are NOT "
160 "considered to be peaks."
161 "This property is enabled when the PeakFindingStrategy has "
162 "been set to NumberOfEventsNormalization.\n"
163 "The value of boxes which contain peaks will be above 1. See "
164 "the below for more information.\n"
168 std::make_unique<EnabledWhenProperty>(
"PeakFindingStrategy",
173 "This will calculate the goniometer rotation (around y-axis "
174 "only) for a constant wavelength. This only works for Q "
175 "sample workspaces.");
177 auto nonNegativeDbl = std::make_shared<BoundedValidator<double>>();
178 nonNegativeDbl->setLower(0);
180 "Wavelength to use when calculating goniometer angle. If not"
181 "set will use the wavelength parameter on the instrument.");
183 std::make_unique<EnabledWhenProperty>(
"CalculateGoniometerForCW",
187 "Whether the goniometer to be calculated is the most inner "
188 "(phi) or most outer (omega)");
190 std::make_unique<EnabledWhenProperty>(
"CalculateGoniometerForCW",
194 "Used when calculating goniometer angle if the q_lab x value "
195 "should be negative, hence the detector of the other side "
196 "(right) of the beam");
199 std::vector<std::string> peakTypes = {
"Automatic",
"Peak",
"LeanElasticPeak"};
200 declareProperty(
"OutputType",
"Automatic", std::make_shared<StringListValidator>(peakTypes),
201 "Type of Peak in OutputWorkspace");
203 "An output PeaksWorkspace with the peaks' found positions.");
206 "If checked, then append the peaks in the output workspace "
208 "If unchecked, the output workspace is replaced (Default).");
210 auto nonNegativeInt = std::make_shared<BoundedValidator<int>>();
211 nonNegativeInt->setLower(0);
212 declareProperty(
"EdgePixels", 0, nonNegativeInt,
"Remove peaks that are at pixels this close to edge. ");
219 inst = ei->getInstrument();
227 }
catch (std::exception &e) {
228 g_log.
warning() <<
"Error finding goniometer matrix. It will not be set in "
229 "the peaks found.\n";
236 std::string dim0 = ws->getDimension(0)->getName();
239 throw std::runtime_error(
"Cannot find peaks in a workspace that is already in HKL space.");
240 }
else if (dim0 ==
"Q_lab_x") {
242 }
else if (dim0 ==
"Q_sample_x")
245 throw std::runtime_error(
"Unexpected dimensions: need either Q_lab_x or Q_sample_x.");
253 if (peakType ==
"Automatic") {
254 if (numExperimentInfo == 0)
256 }
else if (peakType ==
"LeanElasticPeak") {
259 if (numExperimentInfo == 0)
260 throw std::runtime_error(
"Cannot create Peak output with 0 expInfo");
273 auto p = this->
createPeak(Q, binCount, tracer);
278 if (p->getDetectorID() != -1)
280 }
catch (std::exception &e) {
281 g_log.
notice() <<
"Error creating peak at " << Q <<
" because of '" << e.what() <<
"'. Peak will be skipped.\n";
303 std::shared_ptr<DataObjects::Peak> p;
306 p = std::make_shared<Peak>(
inst, Q);
311 bool calcGoniometer =
getProperty(
"CalculateGoniometerForCW");
313 if (calcGoniometer) {
316 if (wavelength == DBL_MAX) {
317 if (
inst->hasParameter(
"wavelength")) {
318 wavelength =
inst->getNumberParameter(
"wavelength").at(0);
320 throw std::runtime_error(
"Could not get wavelength, neither "
321 "Wavelength algorithm property "
322 "set nor instrument wavelength parameter");
329 g_log.
information() <<
"Found goniometer rotation to be in YZY convention [" << angles[0] <<
", " << angles[1]
330 <<
". " << angles[2] <<
"] degrees for Q sample = " << Q <<
"\n";
331 p = std::make_shared<Peak>(
inst, Q, goniometer.
getR());
337 throw std::invalid_argument(
"Cannot Integrate peaks unless the dimension is QLAB or QSAMPLE");
341 p->findDetector(tracer);
345 p->setBinCount(binCount);
354std::shared_ptr<DataObjects::LeanElasticPeak>
356 std::shared_ptr<DataObjects::LeanElasticPeak> p;
361 p = std::make_shared<LeanElasticPeak>(Q);
363 throw std::invalid_argument(
"Cannot find peaks unless the dimension is QSAMPLE");
366 p->setBinCount(binCount);
378 throw std::invalid_argument(
"Workspace must have at least 3 dimensions.");
380 if (isFullMDEvent<MDE, nd>()) {
384 g_log.
warning(
"Workspace contains only lean events. Resultant "
385 "PeaksWorkspaces will not contain full detector "
389 progress(0.01,
"Refreshing Centroids");
401 if (!std::isfinite(threshold)) {
402 g_log.
warning() <<
"Infinite or NaN overall density found. Your input data "
403 "may be invalid. Using a 0 threshold instead.\n";
411 typename std::vector<API::IMDNode *> boxes;
415 ws->
getBox()->getBoxes(boxes, 1000,
true);
418 using dens_box = std::pair<double, API::IMDNode *>;
422 typename std::multimap<double, API::IMDNode *> sortedBoxes;
425 progress(0.20,
"Sorting Boxes by Density");
426 auto it1 = boxes.begin();
427 auto it1_end = boxes.end();
428 for (; it1 != it1_end; it1++) {
433 if (
value > threshold)
434 sortedBoxes.insert(dens_box(
value, box));
439 std::vector<API::IMDNode *> peakBoxes;
444 bool isMDEvent(ws->
id().find(
"MDEventWorkspace") != std::string::npos);
446 int64_t numBoxesFound = 0;
449 typename std::multimap<double, boxPtr>::reverse_iterator it2;
450 auto it2_end = sortedBoxes.rend();
451 for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
453 boxPtr box = it2->second;
454#ifndef MDBOX_TRACK_CENTROID
458 const coord_t *boxCenter = box->getCentroid();
463 for (
auto &peakBoxe : peakBoxes) {
465#ifndef MDBOX_TRACK_CENTROID
467 (*it3)->calculateCentroid(otherCenter);
469 const coord_t *otherCenter = peakBoxe->getCentroid();
474 for (
size_t d = 0;
d < nd;
d++) {
475 coord_t dist = otherCenter[
d] - boxCenter[
d];
476 distSquared += (dist * dist);
489 g_log.
notice() <<
"Number of peaks found exceeded the limit of " <<
m_maxPeaks <<
". Stopping peak finding.\n";
493 peakBoxes.emplace_back(box);
495 for (
size_t d = 0;
d < nd;
d++)
497 g_log.
debug() <<
"; Density = " << density <<
'\n';
499 prog->report(
"Finding Peaks");
503 prog->resetNumSteps(numBoxesFound, 0.95, 1.0);
507 if (numExperimentInfo == 0) {
509 for (
auto box : peakBoxes) {
511#ifndef MDBOX_TRACK_CENTROID
519 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
524 binCount =
static_cast<double>(box->
getNPoints());
530 prog->report(
"Adding Peaks");
540 peakWS->copyExperimentInfoFrom(ei.get());
543 for (
auto box : peakBoxes) {
545 if (numExperimentInfo > 1) {
547 const std::vector<MDE> &events = mdbox->
getEvents();
548 if (std::none_of(events.cbegin(), events.cend(), [&iexp, &numExperimentInfo](
MDE event) {
549 return event.getExpInfoIndex() == iexp || event.getExpInfoIndex() >= numExperimentInfo;
557 if (ei->run().getNumGoniometers() > 1) {
558 const std::vector<MDE> &events =
dynamic_cast<MDBox<MDE, nd> *
>(box)->getEvents();
561 for (
const auto &event : events) {
562 if (event.getExpInfoIndex() == iexp) {
563 sum +=
event.getGoniometerIndex();
571#ifndef MDBOX_TRACK_CENTROID
579 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
585 binCount =
static_cast<double>(box->
getNPoints());
591 auto p = this->
createPeak(Q, binCount, tracer);
595 throw std::runtime_error(
"Failed to cast box to MDBoxBase");
597 addDetectors(*p, *mdBox);
599 if (p->getDetectorID() != -1) {
607 g_log.
information() <<
"Add new peak with Q-center = " << Q[0] <<
", " << Q[1] <<
", " << Q[2] <<
"\n";
609 }
catch (std::exception &e) {
610 g_log.
notice() <<
"Error creating peak at " << Q <<
" because of '" << e.what()
611 <<
"'. Peak will be skipped.\n";
616 prog->report(
"Adding Peaks");
630 size_t nd = ws->getNumDims();
632 throw std::invalid_argument(
"Workspace must have at least 3 dimensions.");
634 g_log.
warning(
"Workspace is an MDHistoWorkspace. Resultant PeaksWorkspaces "
635 "will not contain full detector information.");
638 using dens_box = std::pair<double, size_t>;
642 std::multimap<double, size_t> sortedBoxes;
644 size_t numBoxes = ws->getNPoints();
647 progress(0.10,
"Counting Total Signal");
648 double totalSignal = 0;
649 for (
size_t i = 0; i < numBoxes; i++)
650 totalSignal += ws->getSignalAt(i);
652 double thresholdDensity =
654 if (!std::isfinite(thresholdDensity)) {
655 g_log.
warning() <<
"Infinite or NaN overall density found. Your input data "
656 "may be invalid. Using a 0 threshold instead.\n";
657 thresholdDensity = 0;
659 g_log.
information() <<
"Threshold signal density: " << thresholdDensity <<
'\n';
662 progress(0.20,
"Sorting Boxes by Density");
663 for (
size_t i = 0; i < numBoxes; i++) {
666 if (density > thresholdDensity)
667 sortedBoxes.insert(dens_box(density, i));
672 std::vector<size_t> peakBoxes;
676 int64_t numBoxesFound = 0;
679 std::multimap<double, size_t>::reverse_iterator it2;
680 auto it2_end = sortedBoxes.rend();
681 for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
683 size_t index = it2->second;
685 VMD boxCenter = ws->getCenter(
index);
689 for (
auto &peakBoxe : peakBoxes) {
690 VMD otherCenter = ws->getCenter(peakBoxe);
694 for (
size_t d = 0;
d < nd;
d++) {
695 coord_t dist = otherCenter[
d] - boxCenter[
d];
696 distSquared += (dist * dist);
709 g_log.
notice() <<
"Number of peaks found exceeded the limit of " <<
m_maxPeaks <<
". Stopping peak finding.\n";
713 peakBoxes.emplace_back(
index);
715 g_log.
debug() <<
"; Density = " << density <<
'\n';
717 prog->report(
"Finding Peaks");
721 if (ws->getNumExperimentInfo() == 0) {
723 for (
auto index : peakBoxes) {
725 VMD boxCenter = ws->getCenter(
index);
728 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
737 prog->report(
"Adding Peaks");
741 for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
747 peakWS->copyExperimentInfoFrom(ei.get());
750 for (
auto index : peakBoxes) {
752 VMD boxCenter = ws->getCenter(
index);
755 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
767 prog->report(
"Adding Peaks");
788 uint16_t numExperimentInfo = 0;
790 numExperimentInfo = inMDHW->getNumExperimentInfo();
792 numExperimentInfo = inMDEW->getNumExperimentInfo();
801 if (!
peakWS || !AppendPeaks) {
809 double PeakDistanceThreshold =
getProperty(
"PeakDistanceThreshold");
815 std::string strategy =
getProperty(
"PeakFindingStrategy");
827 throw std::runtime_error(
"This algorithm can only find peaks on a "
828 "MDHistoWorkspace or a MDEventWorkspace; it does "
829 "not work on a regular MatrixWorkspace.");
833 std::vector<std::pair<std::string, bool>> criteria;
834 criteria.emplace_back(
"RunNumber",
true);
835 auto isPeaksWorkspace = std::dynamic_pointer_cast<PeaksWorkspace>(
peakWS);
836 if (isPeaksWorkspace)
837 criteria.emplace_back(
"BankName",
true);
838 criteria.emplace_back(
"bincount",
false);
841 for (
auto i = 0; i !=
peakWS->getNumberPeaks(); ++i) {
853 std::map<std::string, std::string> result;
855 std::string strategy =
getProperty(
"PeakFindingStrategy");
862 if (useNumberOfEventsNormalization && !inMDEW) {
863 result[
"PeakFindingStrategy"] =
"The NumberOfEventsNormalization selection "
864 "can only be used with an MDEventWorkspace "
868 uint16_t numExperimentInfo = 0;
870 numExperimentInfo = inMDHW->getNumExperimentInfo();
872 numExperimentInfo = inMDEW->getNumExperimentInfo();
876 if (peakType ==
"Peak" && numExperimentInfo == 0)
877 result[
"OutputType"] =
"The InputWorkspace doesn't contain any experiment information so the "
878 "OutputType cannot be Peak.";
#define DECLARE_ALGORITHM(classname)
double value
The value of the point.
std::map< DeltaEMode::Type, std::string > index
#define CALL_MDEVENT_FUNCTION3(funcname, workspace)
Macro that makes it possible to call a templated method for a MDEventWorkspace using a IMDEventWorksp...
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.
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
virtual void calculateCentroid(coord_t *) const =0
Calculate the centroid of this box and all sub-boxes.
virtual signal_t getSignalByNEvents() const
virtual size_t getNumChildren() const =0
Get the # of children MDBoxBase'es (non-recursive)
virtual coord_t * getCentroid() const =0
Get the centroid of this box and all sub-boxes.
virtual uint64_t getNPoints() const =0
Get total number of points both in memory and on file if present;.
uint16_t getNumExperimentInfo() const
ExperimentInfo_sptr getExperimentInfo(const uint16_t expInfoIndex)
Get the ExperimentInfo for the given Experiment-Info Index.
A property class for workspaces.
The class LeanElasticPeaksWorkspace stores information about a set of SCD lean peaks.
Templated super-class of a multi-dimensional event "box".
signal_t getSignalNormalized() const override
Returns the integrated signal from all points within, normalized for the cell volume.
Templated class for a multi-dimensional event "box".
std::vector< MDE > & getEvents()
Get vector of events to change.
std::shared_ptr< MDEventWorkspace< MDE, nd > > sptr
Typedef for a shared pointer of this kind of event workspace.
MDBoxBase< MDE, nd > * getBox()
const std::string id() const override
A string ID for the class.
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Structure describing a single-crystal peak.
void addContributingDetID(const int id)
Add a detector ID that contributed to this peak.
The class PeaksWorkspace stores information about a set of SCD peaks.
Class to represent a particular goniometer setting, which is described by the rotation matrix.
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
const Kernel::DblMatrix & getR() const
Return global rotation matrix.
void calcFromQSampleAndWavelength(const Mantid::Kernel::V3D &Q, double wavelength, bool flip_x=false, bool inner=false)
Calculate goniometer for rotation around y-asix for constant wavelength from Q Sample.
Structure describing a single-crystal peak.
virtual void setPeakNumber(int m_PeakNumber)=0
This class is responsible for tracking rays and accumulating a list of objects that are intersected a...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertySettings(const std::string &name, std::unique_ptr< IPropertySettings > settings)
void debug(const std::string &msg)
Logs at debug level.
void notice(const std::string &msg)
Logs at notice level.
void warning(const std::string &msg)
Logs at warning level.
void information(const std::string &msg)
Logs at information level.
The concrete, templated class for properties.
Mantid::Geometry::Instrument_const_sptr inst
Instrument.
coord_t peakRadiusSquared
Estimated radius of peaks. Boxes closer than this are rejected.
void addPeak(const Mantid::Kernel::V3D &Q, const double binCount, const Geometry::InstrumentRayTracer &tracer)
Adds a peak based on Q, bin count & a set of detector IDs.
eDimensionType dimType
Dimension type.
double m_signalThresholdFactor
Signal density factor.
bool m_useNumberOfEventsNormalization
Use number of events normalization for event workspaces.
std::shared_ptr< DataObjects::Peak > createPeak(const Mantid::Kernel::V3D &Q, const double binCount, const Geometry::InstrumentRayTracer &tracer)
Adds a peak based on Q, bin count.
void init() override
Initialise the properties.
void findPeaksHisto(const Mantid::DataObjects::MDHistoWorkspace_sptr &ws)
Run find peaks on a histo workspace.
void readExperimentInfo(const Mantid::API::ExperimentInfo_sptr &ei)
Read member variables from experiment info.
void exec() override
Run the algorithm.
void checkWorkspaceDims(const Mantid::API::IMDWorkspace_sptr &ws)
double DensityThresholdFactor
Thresholding factor.
signal_t m_densityScaleFactor
Arbitrary scaling factor for density to make more manageable numbers, especially for older file forma...
void determineOutputType(const std::string &peakType, const uint16_t numExperimentInfo)
FindPeaksMD()
Constructor.
int64_t m_maxPeaks
Max # of peaks.
void findPeaks(typename DataObjects::MDEventWorkspace< MDE, nd >::sptr ws)
Run find peaks on an MDEventWorkspace.
bool m_addDetectors
Flag to include the detectors within the peak.
std::unique_ptr< Mantid::API::Progress > prog
Progress reporter.
Mantid::API::IPeaksWorkspace_sptr peakWS
Output PeaksWorkspace.
static const std::string numberOfEventsNormalization
NumberOfEventNormalization.
std::map< std::string, std::string > validateInputs() override
Validate the inputs.
int m_edge
Number of edge pixels with no peaks.
Mantid::Kernel::Matrix< double > m_goniometer
Goniometer matrix.
std::shared_ptr< DataObjects::LeanElasticPeak > createLeanElasticPeak(const Mantid::Kernel::V3D &Q, const double binCount, const bool useGoniometer=false)
Adds a peak based on Q, bin count.
int m_runNumber
Run number of the peaks.
static const std::string volumeNormalization
VolumeNormalization.
void addLeanElasticPeak(const Mantid::Kernel::V3D &Q, const double binCount, const bool useGoniometer=false)
Adds a peak based on Q, bin count.
std::shared_ptr< IMDEventWorkspace > IMDEventWorkspace_sptr
Shared pointer to Mantid::API::IMDEventWorkspace.
std::shared_ptr< ExperimentInfo > ExperimentInfo_sptr
Shared pointer to ExperimentInfo.
std::shared_ptr< IMDWorkspace > IMDWorkspace_sptr
Shared pointer to the IMDWorkspace base class.
std::shared_ptr< LeanElasticPeaksWorkspace > LeanElasticPeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< MDHistoWorkspace > MDHistoWorkspace_sptr
A shared pointer to a MDHistoWorkspace.
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
MANTID_GEOMETRY_DLL bool edgePixel(const Geometry::Instrument_const_sptr &inst, const std::string &bankName, int col, int row, int Edge)
Function to find peaks near detector edge.
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
double signal_t
Typedef for the signal recorded in a MDBox, etc.
@ Input
An input workspace.
@ Output
An output workspace.