Mantid
Loading...
Searching...
No Matches
FindPeaksMD.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
8#include "MantidAPI/Run.h"
17#include "MantidKernel/VMD.h"
18
19#include <map>
20#include <vector>
21
22using namespace Mantid::Kernel;
23using namespace Mantid::API;
24using namespace Mantid::DataObjects;
25using namespace Mantid::DataObjects;
26
27namespace Mantid::MDAlgorithms {
28namespace {
29// ---------- Template deduction of the event type
30// --------------------------------
31// See boost::type_traits documentation
32
34template <typename MDE, size_t nd> struct IsFullEvent : boost::false_type {};
35
37template <size_t nd> struct IsFullEvent<MDEvent<nd>, nd> : boost::true_type {};
38
43template <typename MDE, size_t nd> bool isFullMDEvent(const boost::true_type & /*unused*/) { return true; }
44
49template <typename MDE, size_t nd> bool isFullMDEvent(const boost::false_type & /*unused*/) { return false; }
50
54template <typename MDE, size_t nd> bool isFullMDEvent() { return isFullMDEvent<MDE, nd>(IsFullEvent<MDE, nd>()); }
55
61template <typename MDE, size_t nd>
62void addDetectors(DataObjects::Peak &peak, MDBoxBase<MDE, nd> &box, const boost::true_type & /*unused*/) {
63 if (box.getNumChildren() > 0) {
64 std::cerr << "Box has children\n";
65 addDetectors(peak, box, boost::true_type());
66 }
67 auto *mdBox = dynamic_cast<MDBox<MDE, nd> *>(&box);
68 if (!mdBox) {
69 throw std::invalid_argument("FindPeaksMD::addDetectors - Unexpected Box "
70 "type, cannot retrieve events");
71 }
72 const auto &events = mdBox->getConstEvents();
73 auto itend = events.end();
74 for (auto it = events.begin(); it != itend; ++it) {
75 peak.addContributingDetID(it->getDetectorID());
76 }
77}
78
81template <typename MDE, size_t nd>
82void addDetectors(DataObjects::Peak & /*unused*/, MDBoxBase<MDE, nd> & /*unused*/,
83 const boost::false_type & /*unused*/) {
84 throw std::runtime_error("FindPeaksMD - Workspace contains lean events, "
85 "cannot include detector information");
86}
87
93template <typename MDE, size_t nd> void addDetectors(DataObjects::Peak &peak, MDBoxBase<MDE, nd> &box) {
94 // Compile time deduction of the correct function call
95 addDetectors(peak, box, IsFullEvent<MDE, nd>());
96}
97} // namespace
98
99// Register the algorithm into the AlgorithmFactory
100DECLARE_ALGORITHM(FindPeaksMD)
101
102const std::string FindPeaksMD::volumeNormalization = "VolumeNormalization";
103const std::string FindPeaksMD::numberOfEventsNormalization = "NumberOfEventsNormalization";
104
105//----------------------------------------------------------------------------------------------
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() {}
111
112//----------------------------------------------------------------------------------------------
116 declareProperty(std::make_unique<WorkspaceProperty<IMDWorkspace>>("InputWorkspace", "", Direction::Input),
117 "An input MDEventWorkspace or MDHistoWorkspace with at least "
118 "3 dimensions.");
119
120 declareProperty(std::make_unique<PropertyWithValue<double>>("PeakDistanceThreshold", 0.1, Direction::Input),
121 "Threshold distance for rejecting peaks that are found to be too close "
122 "from each other.\n"
123 "This should be some multiple of the radius of a peak. Default: 0.1.");
124
125 declareProperty(std::make_unique<PropertyWithValue<int64_t>>("MaxPeaks", 500, Direction::Input),
126 "Maximum number of peaks to find. Default: 500.");
127
128 std::vector<std::string> strategy = {volumeNormalization, numberOfEventsNormalization};
129 declareProperty("PeakFindingStrategy", volumeNormalization, std::make_shared<StringListValidator>(strategy),
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 "
137 "histogram-based"
138 "raw data which has been converted to an EventWorkspace. The threshold "
139 "for"
140 "peak finding can be controlled by the SingalThresholdFactor property "
141 "which should"
142 "be larger than 1. Note that this approach does not work for event-based "
143 "raw data.\n");
144
145 declareProperty(std::make_unique<PropertyWithValue<double>>("DensityThresholdFactor", 10.0, Direction::Input),
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"
150 "Default: 10.0");
151
152 setPropertySettings("DensityThresholdFactor",
153 std::make_unique<EnabledWhenProperty>(
155
156 declareProperty(std::make_unique<PropertyWithValue<double>>("SignalThresholdFactor", 1.5, Direction::Input),
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"
165 "Default: 1.50");
166
167 setPropertySettings("SignalThresholdFactor",
168 std::make_unique<EnabledWhenProperty>("PeakFindingStrategy",
171
172 declareProperty("CalculateGoniometerForCW", false,
173 "This will calculate the goniometer rotation (around y-axis "
174 "only) for a constant wavelength. This only works for Q "
175 "sample workspaces.");
176
177 auto nonNegativeDbl = std::make_shared<BoundedValidator<double>>();
178 nonNegativeDbl->setLower(0);
179 declareProperty("Wavelength", DBL_MAX, nonNegativeDbl,
180 "Wavelength to use when calculating goniometer angle. If not"
181 "set will use the wavelength parameter on the instrument.");
182 setPropertySettings("Wavelength",
183 std::make_unique<EnabledWhenProperty>("CalculateGoniometerForCW",
185
186 declareProperty("InnerGoniometer", false,
187 "Whether the goniometer to be calculated is the most inner "
188 "(phi) or most outer (omega)");
189 setPropertySettings("InnerGoniometer",
190 std::make_unique<EnabledWhenProperty>("CalculateGoniometerForCW",
192
193 declareProperty("FlipX", false,
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");
197 setPropertySettings("FlipX", std::make_unique<EnabledWhenProperty>(
198 "CalculateGoniometerForCW", Mantid::Kernel::ePropertyCriterion::IS_NOT_DEFAULT));
199 std::vector<std::string> peakTypes = {"Automatic", "Peak", "LeanElasticPeak"};
200 declareProperty("OutputType", "Automatic", std::make_shared<StringListValidator>(peakTypes),
201 "Type of Peak in OutputWorkspace");
202 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("OutputWorkspace", "", Direction::Output),
203 "An output PeaksWorkspace with the peaks' found positions.");
204
205 declareProperty("AppendPeaks", false,
206 "If checked, then append the peaks in the output workspace "
207 "if it exists. \n"
208 "If unchecked, the output workspace is replaced (Default).");
209
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. ");
213}
214
215//----------------------------------------------------------------------------------------------
218 // Instrument associated with workspace
219 inst = ei->getInstrument();
220 // Find the run number
221 m_runNumber = ei->getRunNumber();
222
223 // Find the goniometer rotation matrix
224 m_goniometer = Mantid::Kernel::Matrix<double>(3, 3, true); // Default IDENTITY matrix
225 try {
226 m_goniometer = ei->mutableRun().getGoniometerMatrix();
227 } catch (std::exception &e) {
228 g_log.warning() << "Error finding goniometer matrix. It will not be set in "
229 "the peaks found.\n";
230 g_log.warning() << e.what() << '\n';
231 }
232}
233
235 // Check that the workspace dimensions are in Q-sample-frame or Q-lab-frame.
236 std::string dim0 = ws->getDimension(0)->getName();
237 if (dim0 == "H") {
238 dimType = HKL;
239 throw std::runtime_error("Cannot find peaks in a workspace that is already in HKL space.");
240 } else if (dim0 == "Q_lab_x") {
241 dimType = QLAB;
242 } else if (dim0 == "Q_sample_x")
244 else
245 throw std::runtime_error("Unexpected dimensions: need either Q_lab_x or Q_sample_x.");
246}
247
248void FindPeaksMD::determineOutputType(const std::string &peakType, const uint16_t numExperimentInfo) {
249 // This method will be expanded later to check a property on the
250 // input workspace which can specify a default peak type for that
251 // instrument.
252 m_leanElasticPeak = false;
253 if (peakType == "Automatic") {
254 if (numExperimentInfo == 0)
255 m_leanElasticPeak = true;
256 } else if (peakType == "LeanElasticPeak") {
257 m_leanElasticPeak = true;
258 } else { // Peak
259 if (numExperimentInfo == 0)
260 throw std::runtime_error("Cannot create Peak output with 0 expInfo");
261 }
262}
263
264//----------------------------------------------------------------------------------------------
271void FindPeaksMD::addPeak(const V3D &Q, const double binCount, const Geometry::InstrumentRayTracer &tracer) {
272 try {
273 auto p = this->createPeak(Q, binCount, tracer);
274 if (m_edge > 0) {
275 if (edgePixel(inst, p->getBankName(), p->getCol(), p->getRow(), m_edge))
276 return;
277 }
278 if (p->getDetectorID() != -1)
279 peakWS->addPeak(*p);
280 } catch (std::exception &e) {
281 g_log.notice() << "Error creating peak at " << Q << " because of '" << e.what() << "'. Peak will be skipped.\n";
282 }
283}
284
285//----------------------------------------------------------------------------------------------
293void FindPeaksMD::addLeanElasticPeak(const V3D &Q, const double binCount, const bool useGoniometer) {
294 auto p = this->createLeanElasticPeak(Q, binCount, useGoniometer);
295 peakWS->addPeak(*p);
296}
297
301std::shared_ptr<DataObjects::Peak> FindPeaksMD::createPeak(const Mantid::Kernel::V3D &Q, const double binCount,
302 const Geometry::InstrumentRayTracer &tracer) {
303 std::shared_ptr<DataObjects::Peak> p;
304 if (dimType == QLAB) {
305 // Build using the Q-lab-frame constructor
306 p = std::make_shared<Peak>(inst, Q);
307 // Save gonio matrix for later
308 p->setGoniometerMatrix(m_goniometer);
309 } else if (dimType == QSAMPLE) {
310 // Build using the Q-sample-frame constructor
311 bool calcGoniometer = getProperty("CalculateGoniometerForCW");
312
313 if (calcGoniometer) {
314 // Calculate Q lab from Q sample and wavelength
315 double wavelength = getProperty("Wavelength");
316 if (wavelength == DBL_MAX) {
317 if (inst->hasParameter("wavelength")) {
318 wavelength = inst->getNumberParameter("wavelength").at(0);
319 } else {
320 throw std::runtime_error("Could not get wavelength, neither "
321 "Wavelength algorithm property "
322 "set nor instrument wavelength parameter");
323 }
324 }
325
327 goniometer.calcFromQSampleAndWavelength(Q, wavelength, getProperty("FlipX"), getProperty("InnerGoniometer"));
328 std::vector<double> angles = goniometer.getEulerAngles("YZY");
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());
332
333 } else {
334 p = std::make_shared<Peak>(inst, Q, m_goniometer);
335 }
336 } else {
337 throw std::invalid_argument("Cannot Integrate peaks unless the dimension is QLAB or QSAMPLE");
338 }
339
340 try { // Look for a detector
341 p->findDetector(tracer);
342 } catch (...) { /* Ignore errors in ray-tracer */
343 }
344
345 p->setBinCount(binCount);
346 // Save the run number found before.
347 p->setRunNumber(m_runNumber);
348 return p;
349}
350
354std::shared_ptr<DataObjects::LeanElasticPeak>
355FindPeaksMD::createLeanElasticPeak(const Mantid::Kernel::V3D &Q, const double binCount, const bool useGoniometer) {
356 std::shared_ptr<DataObjects::LeanElasticPeak> p;
357 if (dimType == QSAMPLE) {
358 if (useGoniometer)
359 p = std::make_shared<LeanElasticPeak>(Q, m_goniometer);
360 else
361 p = std::make_shared<LeanElasticPeak>(Q);
362 } else {
363 throw std::invalid_argument("Cannot find peaks unless the dimension is QSAMPLE");
364 }
365
366 p->setBinCount(binCount);
367 // Save the run number found before.
368 p->setRunNumber(m_runNumber);
369 return p;
370}
371//----------------------------------------------------------------------------------------------
376template <typename MDE, size_t nd> void FindPeaksMD::findPeaks(typename MDEventWorkspace<MDE, nd>::sptr ws) {
377 if (nd < 3)
378 throw std::invalid_argument("Workspace must have at least 3 dimensions.");
379
380 if (isFullMDEvent<MDE, nd>()) {
381 m_addDetectors = true;
382 } else {
383 m_addDetectors = false;
384 g_log.warning("Workspace contains only lean events. Resultant "
385 "PeaksWorkspaces will not contain full detector "
386 "information.");
387 }
388
389 progress(0.01, "Refreshing Centroids");
390
391 // TODO: This might be slow, progress report?
392 // Make sure all centroids are fresh
393 // ws->getBox()->refreshCentroid();
394
395 // Calculate a threshold below which a box is too diffuse to be considered a
396 // peak.
397 signal_t threshold = m_useNumberOfEventsNormalization ? ws->getBox()->getSignalByNEvents() * m_signalThresholdFactor
399 threshold *= m_densityScaleFactor;
400
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";
404 threshold = 0;
405 }
406 g_log.information() << "Threshold signal density: " << threshold << '\n';
407
408 using boxPtr = API::IMDNode *;
409 // We will fill this vector with pointers to all the boxes (up to a given
410 // depth)
411 typename std::vector<API::IMDNode *> boxes;
412
413 // Get all the MDboxes
414 progress(0.10, "Getting Boxes");
415 ws->getBox()->getBoxes(boxes, 1000, true);
416
417 // This pair is the <density, ptr to the box>
418 using dens_box = std::pair<double, API::IMDNode *>;
419
420 // Map that will sort the boxes by increasing density. The key = density;
421 // value = box *.
422 typename std::multimap<double, API::IMDNode *> sortedBoxes;
423
424 // --------------- Sort and Filter by Density -----------------------------
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++) {
429 auto box = *it1;
432 // Skip any boxes with too small a signal value.
433 if (value > threshold)
434 sortedBoxes.insert(dens_box(value, box));
435 }
436
437 // --------------- Find Peak Boxes -----------------------------
438 // List of chosen possible peak boxes.
439 std::vector<API::IMDNode *> peakBoxes;
440
441 prog = std::make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
442
443 // used for selecting method for calculating BinCount
444 bool isMDEvent(ws->id().find("MDEventWorkspace") != std::string::npos);
445
446 int64_t numBoxesFound = 0;
447 // Now we go (backwards) through the map
448 // e.g. from highest density down to lowest density.
449 typename std::multimap<double, boxPtr>::reverse_iterator it2;
450 auto it2_end = sortedBoxes.rend();
451 for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
452 signal_t density = it2->first;
453 boxPtr box = it2->second;
454#ifndef MDBOX_TRACK_CENTROID
455 coord_t boxCenter[nd];
456 box->calculateCentroid(boxCenter);
457#else
458 const coord_t *boxCenter = box->getCentroid();
459#endif
460
461 // Compare to all boxes already picked.
462 bool badBox = false;
463 for (auto &peakBoxe : peakBoxes) {
464
465#ifndef MDBOX_TRACK_CENTROID
466 coord_t otherCenter[nd];
467 (*it3)->calculateCentroid(otherCenter);
468#else
469 const coord_t *otherCenter = peakBoxe->getCentroid();
470#endif
471
472 // Distance between this box and a box we already put in.
473 coord_t distSquared = 0.0;
474 for (size_t d = 0; d < nd; d++) {
475 coord_t dist = otherCenter[d] - boxCenter[d];
476 distSquared += (dist * dist);
477 }
478
479 // Reject this box if it is too close to another previously found box.
480 if (distSquared < peakRadiusSquared) {
481 badBox = true;
482 break;
483 }
484 }
485
486 // The box was not rejected for another reason.
487 if (!badBox) {
488 if (numBoxesFound++ >= m_maxPeaks) {
489 g_log.notice() << "Number of peaks found exceeded the limit of " << m_maxPeaks << ". Stopping peak finding.\n";
490 break;
491 }
492
493 peakBoxes.emplace_back(box);
494 g_log.debug() << "Found box at ";
495 for (size_t d = 0; d < nd; d++)
496 g_log.debug() << (d > 0 ? "," : "") << boxCenter[d];
497 g_log.debug() << "; Density = " << density << '\n';
498 // Report progres for each box found.
499 prog->report("Finding Peaks");
500 }
501 }
502
503 prog->resetNumSteps(numBoxesFound, 0.95, 1.0);
504
505 uint16_t numExperimentInfo = ws->getNumExperimentInfo();
506
507 if (numExperimentInfo == 0) {
508 // --- Convert the "boxes" to peaks ----
509 for (auto box : peakBoxes) {
510 // The center of the box = Q in the lab frame
511#ifndef MDBOX_TRACK_CENTROID
512 coord_t boxCenter[nd];
513 box->calculateCentroid(boxCenter);
514#else
515 const coord_t *boxCenter = box->getCentroid();
516#endif
517
518 // Q of the centroid of the box
519 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
520
521 // The "bin count" used will be the box density.
522 double binCount = box->getSignalNormalized() * m_densityScaleFactor;
523 if (isMDEvent)
524 binCount = static_cast<double>(box->getNPoints());
525
526 // Create the peak
527 addLeanElasticPeak(Q, binCount);
528
529 // Report progres for each box found.
530 prog->report("Adding Peaks");
531
532 } // for each box found
533 } else {
534 for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
536 this->readExperimentInfo(ei);
537
539 // Copy the instrument, sample, run to the peaks workspace.
540 peakWS->copyExperimentInfoFrom(ei.get());
541
542 // --- Convert the "boxes" to peaks ----
543 for (auto box : peakBoxes) {
544 // If no events from this experimental contribute to the box then skip
545 if (numExperimentInfo > 1) {
546 auto *mdbox = dynamic_cast<MDBox<MDE, nd> *>(box);
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;
550 }))
551 continue;
552 }
553
554 // If multiple goniometers than use the average one from the
555 // events in the box, that matches this expInfoIndex, this assumes
556 // the events are added in same order as the goniometers
557 if (ei->run().getNumGoniometers() > 1) {
558 const std::vector<MDE> &events = dynamic_cast<MDBox<MDE, nd> *>(box)->getEvents();
559 double sum = 0;
560 double count = 0;
561 for (const auto &event : events) {
562 if (event.getExpInfoIndex() == iexp) {
563 sum += event.getGoniometerIndex();
564 count++;
565 }
566 }
567 m_goniometer = ei->mutableRun().getGoniometerMatrix(lrint(sum / count));
568 }
569 // The center of the box = Q in the lab frame
570
571#ifndef MDBOX_TRACK_CENTROID
572 coord_t boxCenter[nd];
573 box->calculateCentroid(boxCenter);
574#else
575 const coord_t *boxCenter = box->getCentroid();
576#endif
577
578 // Q of the centroid of the box
579 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
580
581 // The "bin count" used will be the box density or the number of events
582 // in the box
583 double binCount = box->getSignalNormalized() * m_densityScaleFactor;
584 if (isMDEvent)
585 binCount = static_cast<double>(box->getNPoints());
586
587 if (m_leanElasticPeak) {
588 addLeanElasticPeak(Q, binCount, true);
589 } else {
590 try {
591 auto p = this->createPeak(Q, binCount, tracer);
592 if (m_addDetectors) {
593 auto mdBox = dynamic_cast<MDBoxBase<MDE, nd> *>(box);
594 if (!mdBox) {
595 throw std::runtime_error("Failed to cast box to MDBoxBase");
596 }
597 addDetectors(*p, *mdBox);
598 }
599 if (p->getDetectorID() != -1) {
600 if (m_edge > 0) {
601 if (!edgePixel(inst, p->getBankName(), p->getCol(), p->getRow(), m_edge))
602 peakWS->addPeak(*p);
603 ;
604 } else {
605 peakWS->addPeak(*p);
606 }
607 g_log.information() << "Add new peak with Q-center = " << Q[0] << ", " << Q[1] << ", " << Q[2] << "\n";
608 }
609 } catch (std::exception &e) {
610 g_log.notice() << "Error creating peak at " << Q << " because of '" << e.what()
611 << "'. Peak will be skipped.\n";
612 }
613 }
614
615 // Report progress for each box found.
616 prog->report("Adding Peaks");
617
618 } // for each box found
619 }
620 }
621 g_log.notice() << "Number of peaks found: " << peakWS->getNumberPeaks() << '\n';
622}
623
624//----------------------------------------------------------------------------------------------
630 size_t nd = ws->getNumDims();
631 if (nd < 3)
632 throw std::invalid_argument("Workspace must have at least 3 dimensions.");
633
634 g_log.warning("Workspace is an MDHistoWorkspace. Resultant PeaksWorkspaces "
635 "will not contain full detector information.");
636
637 // This pair is the <density, box index>
638 using dens_box = std::pair<double, size_t>;
639
640 // Map that will sort the boxes by increasing density. The key = density;
641 // value = box index.
642 std::multimap<double, size_t> sortedBoxes;
643
644 size_t numBoxes = ws->getNPoints();
645
646 // --------- Count the overall signal density -----------------------------
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);
651 // Calculate the threshold density
652 double thresholdDensity =
653 (totalSignal * ws->getInverseVolume() / double(numBoxes)) * DensityThresholdFactor * m_densityScaleFactor;
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;
658 }
659 g_log.information() << "Threshold signal density: " << thresholdDensity << '\n';
660
661 // -------------- Sort and Filter by Density -----------------------------
662 progress(0.20, "Sorting Boxes by Density");
663 for (size_t i = 0; i < numBoxes; i++) {
664 double density = ws->getSignalNormalizedAt(i) * m_densityScaleFactor;
665 // Skip any boxes with too small a signal density.
666 if (density > thresholdDensity)
667 sortedBoxes.insert(dens_box(density, i));
668 }
669
670 // --------------- Find Peak Boxes -----------------------------
671 // List of chosen possible peak boxes.
672 std::vector<size_t> peakBoxes;
673
674 prog = std::make_unique<Progress>(this, 0.30, 0.95, m_maxPeaks);
675
676 int64_t numBoxesFound = 0;
677 // Now we go (backwards) through the map
678 // e.g. from highest density down to lowest density.
679 std::multimap<double, size_t>::reverse_iterator it2;
680 auto it2_end = sortedBoxes.rend();
681 for (it2 = sortedBoxes.rbegin(); it2 != it2_end; ++it2) {
682 signal_t density = it2->first;
683 size_t index = it2->second;
684 // Get the center of the box
685 VMD boxCenter = ws->getCenter(index);
686
687 // Compare to all boxes already picked.
688 bool badBox = false;
689 for (auto &peakBoxe : peakBoxes) {
690 VMD otherCenter = ws->getCenter(peakBoxe);
691
692 // Distance between this box and a box we already put in.
693 coord_t distSquared = 0.0;
694 for (size_t d = 0; d < nd; d++) {
695 coord_t dist = otherCenter[d] - boxCenter[d];
696 distSquared += (dist * dist);
697 }
698
699 // Reject this box if it is too close to another previously found box.
700 if (distSquared < peakRadiusSquared) {
701 badBox = true;
702 break;
703 }
704 }
705
706 // The box was not rejected for another reason.
707 if (!badBox) {
708 if (numBoxesFound++ >= m_maxPeaks) {
709 g_log.notice() << "Number of peaks found exceeded the limit of " << m_maxPeaks << ". Stopping peak finding.\n";
710 break;
711 }
712
713 peakBoxes.emplace_back(index);
714 g_log.debug() << "Found box at index " << index;
715 g_log.debug() << "; Density = " << density << '\n';
716 // Report progres for each box found.
717 prog->report("Finding Peaks");
718 }
719 }
720
721 if (ws->getNumExperimentInfo() == 0) {
722 // --- Convert the "boxes" to peaks ----
723 for (auto index : peakBoxes) {
724 // The center of the box = Q in the lab frame
725 VMD boxCenter = ws->getCenter(index);
726
727 // Q of the centroid of the box
728 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
729
730 // The "bin count" used will be the box density.
731 double binCount = ws->getSignalNormalizedAt(index) * m_densityScaleFactor;
732
733 // Create the peak
734 addLeanElasticPeak(Q, binCount);
735
736 // Report progres for each box found.
737 prog->report("Adding Peaks");
738
739 } // for each box found
740 } else {
741 for (uint16_t iexp = 0; iexp < ws->getNumExperimentInfo(); iexp++) {
742 ExperimentInfo_sptr ei = ws->getExperimentInfo(iexp);
743 this->readExperimentInfo(ei);
745
746 // Copy the instrument, sample, run to the peaks workspace.
747 peakWS->copyExperimentInfoFrom(ei.get());
748
749 // --- Convert the "boxes" to peaks ----
750 for (auto index : peakBoxes) {
751 // The center of the box = Q in the lab frame
752 VMD boxCenter = ws->getCenter(index);
753
754 // Q of the centroid of the box
755 V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);
756
757 // The "bin count" used will be the box density.
758 double binCount = ws->getSignalNormalizedAt(index) * m_densityScaleFactor;
759
760 // Create the peak
762 addLeanElasticPeak(Q, binCount, true);
763 else
764 addPeak(Q, binCount, tracer);
765
766 // Report progres for each box found.
767 prog->report("Adding Peaks");
768
769 } // for each box found
770 }
771 }
772 g_log.notice() << "Number of peaks found: " << peakWS->getNumberPeaks() << '\n';
773} // namespace MDAlgorithms
774
775//----------------------------------------------------------------------------------------------
779
780 // The MDEventWorkspace as input
781 IMDWorkspace_sptr inWS = getProperty("InputWorkspace");
782 checkWorkspaceDims(inWS);
783 MDHistoWorkspace_sptr inMDHW = std::dynamic_pointer_cast<MDHistoWorkspace>(inWS);
784 IMDEventWorkspace_sptr inMDEW = std::dynamic_pointer_cast<IMDEventWorkspace>(inWS);
785
786 bool AppendPeaks = getProperty("AppendPeaks");
787
788 uint16_t numExperimentInfo = 0;
789 if (inMDHW)
790 numExperimentInfo = inMDHW->getNumExperimentInfo();
791 else if (inMDEW)
792 numExperimentInfo = inMDEW->getNumExperimentInfo();
793
794 std::string peakType = getProperty("OutputType");
795
796 determineOutputType(peakType, numExperimentInfo);
797
798 // Output peaks workspace, create if needed
799 peakWS = getProperty("OutputWorkspace");
800
801 if (!peakWS || !AppendPeaks) {
804 else
806 }
807
808 // Other parameters
809 double PeakDistanceThreshold = getProperty("PeakDistanceThreshold");
810 peakRadiusSquared = static_cast<coord_t>(PeakDistanceThreshold * PeakDistanceThreshold);
811
812 DensityThresholdFactor = getProperty("DensityThresholdFactor");
813 m_signalThresholdFactor = getProperty("SignalThresholdFactor");
814
815 std::string strategy = getProperty("PeakFindingStrategy");
817
818 m_maxPeaks = getProperty("MaxPeaks");
819 m_edge = this->getProperty("EdgePixels");
820
821 // Execute the proper algo based on the type of workspace
822 if (inMDHW) {
823 this->findPeaksHisto(inMDHW);
824 } else if (inMDEW) {
825 CALL_MDEVENT_FUNCTION3(this->findPeaks, inMDEW);
826 } else {
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.");
830 }
831
832 // Do a sort by bank name (if peaks of class "Peak") and then descending bin count (intensity)
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);
839 peakWS->sort(criteria);
840
841 for (auto i = 0; i != peakWS->getNumberPeaks(); ++i) {
842 Mantid::Geometry::IPeak &p = peakWS->getPeak(i);
843 p.setPeakNumber(i + 1);
844 }
845 // Save the output
846 setProperty("OutputWorkspace", peakWS);
847}
848
849//----------------------------------------------------------------------------------------------
852std::map<std::string, std::string> FindPeaksMD::validateInputs() {
853 std::map<std::string, std::string> result;
854 // Check for number of event normalzation
855 std::string strategy = getProperty("PeakFindingStrategy");
856
857 const bool useNumberOfEventsNormalization = strategy == numberOfEventsNormalization;
858 IMDWorkspace_sptr inWS = getProperty("InputWorkspace");
859 IMDEventWorkspace_sptr inMDEW = std::dynamic_pointer_cast<IMDEventWorkspace>(inWS);
860 MDHistoWorkspace_sptr inMDHW = std::dynamic_pointer_cast<MDHistoWorkspace>(inWS);
861
862 if (useNumberOfEventsNormalization && !inMDEW) {
863 result["PeakFindingStrategy"] = "The NumberOfEventsNormalization selection "
864 "can only be used with an MDEventWorkspace "
865 "as the input.";
866 }
867
868 uint16_t numExperimentInfo = 0;
869 if (inMDHW)
870 numExperimentInfo = inMDHW->getNumExperimentInfo();
871 else if (inMDEW)
872 numExperimentInfo = inMDEW->getNumExperimentInfo();
873
874 std::string peakType = getProperty("OutputType");
875
876 if (peakType == "Peak" && numExperimentInfo == 0)
877 result["OutputType"] = "The InputWorkspace doesn't contain any experiment information so the "
878 "OutputType cannot be Peak.";
879
880 return result;
881}
882
883} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
std::map< DeltaEMode::Type, std::string > index
Definition: DeltaEMode.cpp:19
#define CALL_MDEVENT_FUNCTION3(funcname, workspace)
Macro that makes it possible to call a templated method for a MDEventWorkspace using a IMDEventWorksp...
int count
counter
Definition: Matrix.cpp:37
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
Definition: Algorithm.cpp:1913
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
Kernel::Logger & g_log
Definition: Algorithm.h:451
void progress(double p, const std::string &msg="", double estimatedTime=0.0, int progressPrecision=0)
Sends ProgressNotification.
Definition: Algorithm.cpp:231
virtual void calculateCentroid(coord_t *) const =0
Calculate the centroid of this box and all sub-boxes.
virtual signal_t getSignalByNEvents() const
Definition: IMDNode.h:266
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;.
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".
Definition: MDBoxBase.h:50
signal_t getSignalNormalized() const override
Returns the integrated signal from all points within, normalized for the cell volume.
Definition: MDBoxBase.h:274
Templated class for a multi-dimensional event "box".
Definition: MDBox.h:45
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.
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,...
Definition: MDEvent.h:36
Templated class holding data about a neutron detection event in N-dimensions (for example,...
Definition: MDLeanEvent.h:60
Structure describing a single-crystal peak.
Definition: Peak.h:34
void addContributingDetID(const int id)
Add a detector ID that contributed to this peak.
Definition: Peak.cpp:278
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.
Definition: Goniometer.h:55
std::vector< double > getEulerAngles(const std::string &convention="YZX")
Return Euler angles acording to a convention.
Definition: Goniometer.cpp:282
const Kernel::DblMatrix & getR() const
Return global rotation matrix.
Definition: Goniometer.cpp:80
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.
Definition: Goniometer.cpp:197
Structure describing a single-crystal peak.
Definition: IPeak.h:26
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.
Definition: Logger.cpp:114
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Numerical Matrix class.
Definition: Matrix.h:42
The concrete, templated class for properties.
Class for 3D vectors.
Definition: V3D.h:34
Mantid::Geometry::Instrument_const_sptr inst
Instrument.
Definition: FindPeaksMD.h:111
coord_t peakRadiusSquared
Estimated radius of peaks. Boxes closer than this are rejected.
Definition: FindPeaksMD.h:86
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.
Definition: FindPeaksMD.h:115
double m_signalThresholdFactor
Signal density factor.
Definition: FindPeaksMD.h:122
bool m_useNumberOfEventsNormalization
Use number of events normalization for event workspaces.
Definition: FindPeaksMD.h:120
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.
Definition: FindPeaksMD.h:89
signal_t m_densityScaleFactor
Arbitrary scaling factor for density to make more manageable numbers, especially for older file forma...
Definition: FindPeaksMD.h:102
void determineOutputType(const std::string &peakType, const uint16_t numExperimentInfo)
int64_t m_maxPeaks
Max # of peaks.
Definition: FindPeaksMD.h:92
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.
Definition: FindPeaksMD.h:98
std::unique_ptr< Mantid::API::Progress > prog
Progress reporter.
Definition: FindPeaksMD.h:105
Mantid::API::IPeaksWorkspace_sptr peakWS
Output PeaksWorkspace.
Definition: FindPeaksMD.h:83
static const std::string numberOfEventsNormalization
NumberOfEventNormalization.
Definition: FindPeaksMD.h:126
std::map< std::string, std::string > validateInputs() override
Validate the inputs.
int m_edge
Number of edge pixels with no peaks.
Definition: FindPeaksMD.h:95
Mantid::Kernel::Matrix< double > m_goniometer
Goniometer matrix.
Definition: FindPeaksMD.h:117
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.
Definition: FindPeaksMD.h:113
static const std::string volumeNormalization
VolumeNormalization.
Definition: FindPeaksMD.h:124
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.
Definition: IMDWorkspace.h:146
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.
Definition: EdgePixel.cpp:22
float coord_t
Typedef for the data type to use for coordinate axes in MD objects such as MDBox, MDEventWorkspace,...
Definition: MDTypes.h:27
double signal_t
Typedef for the signal recorded in a MDBox, etc.
Definition: MDTypes.h:36
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54