Mantid
Loading...
Searching...
No Matches
IntegrateEllipsoidsV2.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 +
11#include "MantidAPI/Run.h"
12#include "MantidAPI/Sample.h"
23#include "MantidHistogramData/LinearGenerator.h"
32
33#include <cmath>
34
35using namespace Mantid::API;
36using namespace Mantid::HistogramData;
37using namespace Mantid::Kernel;
38using namespace Mantid::Geometry;
39using namespace Mantid::DataObjects;
40
41namespace Mantid::MDAlgorithms {
42
43// Register the algorithm into the AlgorithmFactory
44DECLARE_ALGORITHM(IntegrateEllipsoidsV2)
45
46
47const std::string ELASTIC("Elastic");
48
50const std::string Q3D("Q3D");
51
53const std::size_t DIMS(3);
54
56 EventWorkspace_sptr &wksp) {
57 auto numSpectra = static_cast<int>(wksp->getNumberHistograms());
59 for (int i = 0; i < numSpectra; ++i) {
61
62 // units conversion helper
63 UnitsConversionHelper unitConverter;
64 unitConverter.initialize(m_targWSDescr, "Momentum");
65
66 // initialize the MD coordinates conversion class
67 MDTransfQ3D qConverter;
68 qConverter.initialize(m_targWSDescr);
69
70 std::vector<double> buffer(DIMS);
71 // get a reference to the event list
72 EventList &events = wksp->getSpectrum(i);
73
75 events.compressEvents(1e-5, &events);
76
77 // check to see if the event list is empty
78 if (events.empty()) {
79 prog.report();
80 continue; // nothing to do
81 }
82
83 // update which pixel is being converted
84 std::vector<Mantid::coord_t> locCoord(DIMS, 0.);
85 unitConverter.updateConversion(i);
86 qConverter.calcYDepCoordinates(locCoord, i);
87
88 // loop over the events
89 double signal(1.); // ignorable garbage
90 double errorSq(1.); // ignorable garbage
91 const std::vector<WeightedEventNoTime> &raw_events = events.getWeightedEventsNoTime();
92 std::vector<std::pair<std::pair<double, double>, V3D>> qList;
93 for (const auto &raw_event : raw_events) {
94 double val = unitConverter.convertUnits(raw_event.tof());
95 qConverter.calcMatrixCoord(val, locCoord, signal, errorSq);
96 for (size_t dim = 0; dim < DIMS; ++dim) {
97 buffer[dim] = locCoord[dim];
98 }
99 V3D qVec(buffer[0], buffer[1], buffer[2]);
100 qList.emplace_back(std::pair<double, double>(raw_event.m_weight, raw_event.m_errorSquared), qVec);
101 } // end of loop over events in list
102 PARALLEL_CRITICAL(addEvents) { integrator.addEvents(qList); }
103
104 prog.report();
106 } // end of loop over spectra
108 integrator.populateCellsWithPeaks();
109}
110
112 auto numSpectra = static_cast<int>(wksp->getNumberHistograms());
114 for (int i = 0; i < numSpectra; ++i) {
116
117 // units conversion helper
118 UnitsConversionHelper unitConverter;
119 unitConverter.initialize(m_targWSDescr, "Momentum");
120
121 // initialize the MD coordinates conversion class
122 MDTransfQ3D qConverter;
123 qConverter.initialize(m_targWSDescr);
124
125 std::vector<double> buffer(DIMS);
126 // get tof and y values
127 const auto &xVals = wksp->points(i);
128 const auto &yVals = wksp->y(i);
129 const auto &eVals = wksp->e(i);
130
131 // update which pixel is being converted
132 std::vector<Mantid::coord_t> locCoord(DIMS, 0.);
133 unitConverter.updateConversion(i);
134 qConverter.calcYDepCoordinates(locCoord, i);
135
136 // loop over the events
137 double signal(1.); // ignorable garbage
138 double errorSq(1.); // ignorable garbage
139
140 SlimEvents qList;
141 for (size_t j = 0; j < yVals.size(); ++j) {
142 const double &yVal = yVals[j];
143 const double &esqVal = eVals[j] * eVals[j]; // error squared (variance)
144 if (yVal > 0) // TODO, is this condition right?
145 {
146 double val = unitConverter.convertUnits(xVals[j]);
147 qConverter.calcMatrixCoord(val, locCoord, signal, errorSq);
148 for (size_t dim = 0; dim < DIMS; ++dim) {
149 buffer[dim] = locCoord[dim]; // TODO. Looks un-necessary to me. Can't
150 // we just add localCoord directly to
151 // qVec
152 }
153 V3D qVec(buffer[0], buffer[1], buffer[2]);
154 if (std::isnan(qVec[0]) || std::isnan(qVec[1]) || std::isnan(qVec[2]))
155 continue;
156 // Account for counts in histograms by increasing the qList with the
157 // same q-point
158 qList.emplace_back(std::pair<double, double>(yVal, esqVal), qVec);
159 }
160 }
161 PARALLEL_CRITICAL(addHisto) { integrator.addEvents(qList); }
162 prog.report();
164 } // end of loop over spectra
166 integrator.populateCellsWithPeaks();
167}
168
171 // not used in this version
172 removeProperty("IntegrateInHKL", false);
173 removeProperty("GetUBFromPeaksWorkspace", false);
174
175 initInstance(*this);
176}
178
179 std::shared_ptr<BoundedValidator<double>> mustBePositive(new BoundedValidator<double>());
180 mustBePositive->setLower(0.0);
181
182 // Different defaults
183 alg.removeProperty("SatelliteRegionRadius", false);
184 alg.removeProperty("SatellitePeakSize", false);
185 alg.removeProperty("SatelliteBackgroundInnerSize", false);
186 alg.removeProperty("SatelliteBackgroundOuterSize", false);
187
188 // satellite realted properties
189 alg.declareProperty("SatelliteRegionRadius", EMPTY_DBL(), mustBePositive,
190 "Only events at most this distance from a satellite peak will be considered when integration");
191 alg.declareProperty("SatellitePeakSize", EMPTY_DBL(), mustBePositive,
192 "Half-length of major axis for satellite peak ellipsoid");
193 alg.declareProperty("ShareBackground", false, "Whether to use the same peak background region for satellite peaks.");
194 alg.declareProperty(
195 "SatelliteBackgroundInnerSize", EMPTY_DBL(), mustBePositive,
196 "Half-length of major axis for the inner ellipsoidal surface of background region of the satellite peak");
197 alg.declareProperty(
198 "SatelliteBackgroundOuterSize", EMPTY_DBL(), mustBePositive,
199 "Half-length of major axis for the outer ellipsoidal surface of background region of the satellite peak");
200}
201
207std::map<std::string, std::string> IntegrateEllipsoidsV2::validateInputs() {
208 std::map<std::string, std::string> issues;
209
210 // case 1: specified peak and background must be realisitc
211 double radius_m = getProperty("RegionRadius");
212 bool specify_size = getProperty("SpecifySize");
213 double peak_radius = getProperty("PeakSize");
214 double back_inner_radius = getProperty("BackgroundInnerSize");
215 double back_outer_radius = getProperty("BackgroundOuterSize");
216 if (specify_size) {
217 if (back_outer_radius > radius_m) {
218 issues["SpecifySize"] = "BackgroundOuterSize must be less than or equal to the RegionRadius";
219 }
220 if (back_inner_radius >= back_outer_radius) {
221 issues["SpecifySize"] = "BackgroundInnerSize must be less than BackgroundOuterSize";
222 }
223 if (peak_radius > back_inner_radius) {
224 issues["SpecifySize"] = "PeakSize must be less than or equal to the BackgroundInnerSize";
225 }
226 }
227
228 // case 2: specified satellite peak and background must be realisitc
229 double satellite_radius = (getPointerToProperty("SatelliteRegionRadius")->isDefault())
230 ? getProperty("RegionRadius")
231 : getProperty("SatelliteRegionRadius");
232 double satellite_peak_radius = (getPointerToProperty("SatellitePeakSize")->isDefault())
233 ? getProperty("PeakSize")
234 : getProperty("SatellitePeakSize");
235 double satellite_back_inner_radius = (getPointerToProperty("SatelliteBackgroundInnerSize")->isDefault())
236 ? getProperty("BackgroundInnerSize")
237 : getProperty("SatelliteBackgroundInnerSize");
238 double satellite_back_outer_radius = (getPointerToProperty("SatelliteBackgroundOuterSize")->isDefault())
239 ? getProperty("BackgroundOuterSize")
240 : getProperty("SatelliteBackgroundOuterSize");
241 if (specify_size) {
242 if (satellite_back_outer_radius > satellite_radius) {
243 issues["SpecifySize"] = "SatelliteBackgroundOuterSize must be less than or equal to the SatelliteRegionRadius";
244 }
245 if (satellite_back_inner_radius > satellite_back_outer_radius) {
246 issues["SpecifySize"] = "SatelliteBackgroundInnerSize must be less than SatelliteBackgroundOuterSize";
247 }
248 if (satellite_peak_radius > satellite_back_inner_radius) {
249 issues["SpecifySize"] = "SatellitePeakSize must be less than or equal to the SatelliteBackgroundInnerSize";
250 }
251 }
252
253 // case 3: anything else?
254
255 return issues;
256}
257
259 // get the input workspace
260 MatrixWorkspace_sptr wksp = getProperty("InputWorkspace");
261
262 EventWorkspace_sptr eventWS = std::dynamic_pointer_cast<EventWorkspace>(wksp);
263 Workspace2D_sptr histoWS = std::dynamic_pointer_cast<Workspace2D>(wksp);
264 if (!eventWS && !histoWS) {
265 throw std::runtime_error("IntegrateEllipsoidsV2 needs either a "
266 "EventWorkspace or Workspace2D as input.");
267 }
268
269 // error out if there are not events
270 if (eventWS && eventWS->getNumberEvents() <= 0) {
271 throw std::runtime_error("IntegrateEllipsoidsV2 does not work for empty event lists");
272 }
273
274 PeaksWorkspace_sptr in_peak_ws = getProperty("PeaksWorkspace");
275 if (!in_peak_ws) {
276 throw std::runtime_error("Could not read the peaks workspace");
277 }
278
279 double radius_m = getProperty("RegionRadius");
280 int numSigmas = getProperty("NumSigmas");
281 double cutoffIsigI = getProperty("CutoffIsigI");
282 bool specify_size = getProperty("SpecifySize");
283 double peak_radius = getProperty("PeakSize");
284 double back_inner_radius = getProperty("BackgroundInnerSize");
285 double back_outer_radius = getProperty("BackgroundOuterSize");
286 bool integrateEdge = getProperty("IntegrateIfOnEdge");
287 bool adaptiveQBackground = getProperty("AdaptiveQBackground");
288 double adaptiveQMultiplier = getProperty("AdaptiveQMultiplier");
289 double adaptiveQBackgroundMultiplier = 0.0;
290 bool useOnePercentBackgroundCorrection = getProperty("UseOnePercentBackgroundCorrection");
291 // satellite related properties
292 // NOTE: fallback to Brag Peak properties if satellite peak related properties are not specified
293 double satellite_radius = (getPointerToProperty("SatelliteRegionRadius")->isDefault())
294 ? getProperty("RegionRadius")
295 : getProperty("SatelliteRegionRadius");
296 double satellite_peak_radius = (getPointerToProperty("SatellitePeakSize")->isDefault())
297 ? getProperty("PeakSize")
298 : getProperty("SatellitePeakSize");
299 double satellite_back_inner_radius = (getPointerToProperty("SatelliteBackgroundInnerSize")->isDefault())
300 ? getProperty("BackgroundInnerSize")
301 : getProperty("SatelliteBackgroundInnerSize");
302 double satellite_back_outer_radius = (getPointerToProperty("SatelliteBackgroundOuterSize")->isDefault())
303 ? getProperty("BackgroundOuterSize")
304 : getProperty("SatelliteBackgroundOuterSize");
305 bool shareBackground = getProperty("ShareBackground");
306
307 if (adaptiveQBackground)
308 adaptiveQBackgroundMultiplier = adaptiveQMultiplier;
309 if (!integrateEdge) {
310 // This only fails in the unit tests which say that MaskBTP is not
311 // registered
312 try {
313 runMaskDetectors(in_peak_ws, "Tube", "edges");
314 runMaskDetectors(in_peak_ws, "Pixel", "edges");
315 } catch (...) {
316 g_log.error("Can't execute MaskBTP algorithm for this instrument to set "
317 "edge for IntegrateIfOnEdge option");
318 }
319 calculateE1(in_peak_ws->detectorInfo()); // fill E1Vec for use in detectorQ
320 }
321
322 Mantid::DataObjects::PeaksWorkspace_sptr peak_ws = getProperty("OutputWorkspace");
323 if (peak_ws != in_peak_ws)
324 peak_ws = in_peak_ws->clone();
325
326 // get the list of peak Q's for the integrator
327 std::vector<Peak> &peaks = peak_ws->getPeaks();
328 size_t n_peaks = peak_ws->getNumberPeaks();
329 SlimEvents qList;
330
331 // Note: we skip un-indexed peaks if index count greater than zero
332 for (size_t i = 0; i < n_peaks; i++) {
333 // check if peak is satellite peak
334 const bool isSatellitePeak = (peaks[i].getIntMNP().norm2() > 0);
335 const V3D peak_q = peaks[i].getQLabFrame();
336 const bool isOrigin = isSatellitePeak ? IntegrateQLabEvents::isOrigin(peak_q, satellite_radius)
337 : IntegrateQLabEvents::isOrigin(peak_q, radius_m);
338 if (isOrigin) {
339 continue; // skip this peak
340 }
341 // add peak Q to list
342 qList.emplace_back(std::pair<double, double>(1., 1.), peak_q);
343 }
344
345 // Peak vectors
346 std::vector<double> PeakRadiusVector(n_peaks, peak_radius);
347 std::vector<double> BackgroundInnerRadiusVector(n_peaks, back_inner_radius);
348 std::vector<double> BackgroundOuterRadiusVector(n_peaks, back_outer_radius);
349 // Satellite peak vectors
350 std::vector<double> SatellitePeakRadiusVector(n_peaks, satellite_peak_radius);
351 std::vector<double> SatelliteBackgroundInnerRadiusVector(n_peaks, satellite_back_inner_radius);
352 std::vector<double> SatelliteBackgroundOuterRadiusVector(n_peaks, satellite_back_outer_radius);
353
354 // make the integrator
355 m_braggPeakRadius = radius_m;
356 m_satellitePeakRadius = satellite_radius;
357
358 IntegrateQLabEvents integrator(qList, satellite_radius, useOnePercentBackgroundCorrection);
359
360 // get the events and add
361 // them to the inegrator
362 // set up a descripter of where we are going
363 this->initTargetWSDescr(wksp);
364
365 // set up the progress bar
366 const size_t numSpectra = wksp->getNumberHistograms();
367 Progress prog(this, 0.5, 1.0, numSpectra);
368
369 // TODO Refactor - Skip this block to find out how many tests will be broken
370 if (eventWS) {
371 // process as EventWorkspace
372 qListFromEventWS(integrator, prog, eventWS);
373 } else {
374 // process as Workspace2D
375 qListFromHistoWS(integrator, prog, histoWS);
376 }
377
378 // map of satellite peaks for each bragg peak
379 std::map<size_t, std::vector<Peak *>> satellitePeakMap;
380 // lists containing indices of bragg or satellite peaks
381 std::vector<size_t> satellitePeaks;
382 if (shareBackground) {
383 pairBraggSatellitePeaks(n_peaks, peaks, satellitePeakMap, satellitePeaks);
384 }
385
386 // Integrate peaks
387 std::vector<double> principalaxis1, principalaxis2, principalaxis3;
388 // cached background and sigma background for each bragg peak (including ellipsoid ratio factor)
389 std::map<size_t, std::pair<double, double>> cachedBraggBackground;
390 for (size_t i = 0; i < n_peaks; i++) {
391 // check if peak is satellite peak
392 const bool isSatellitePeak = (peaks[i].getIntMNP().norm2() > 0);
393 // grab QLabFrame
394 const V3D peak_q = peaks[i].getQLabFrame();
395
396 // check if peak is origin (skip if true)
397 const bool isOrigin = isSatellitePeak ? IntegrateQLabEvents::isOrigin(peak_q, m_satellitePeakRadius)
399 if (isOrigin) {
400 continue;
401 }
402
403 // modulus of Q
404 const double lenQpeak = (adaptiveQMultiplier != 0.0) ? peak_q.norm() : 0.0;
405 // compuate adaptive radius
406 double adaptiveRadius = isSatellitePeak ? adaptiveQMultiplier * lenQpeak + satellite_peak_radius
407 : adaptiveQMultiplier * lenQpeak + peak_radius;
408 // - error checking for adaptive radius
409 if (adaptiveRadius < 0.0) {
410 // Unphysical case: radius is negative
411 std::ostringstream errmsg;
412 errmsg << "Error: Radius for integration sphere of peak " << i << " is negative = " << adaptiveRadius << '\n';
413 g_log.error() << errmsg.str();
414 // zero the peak
415 peaks[i].setIntensity(0.0);
416 peaks[i].setSigmaIntensity(0.0);
417 PeakRadiusVector[i] = 0.0;
418 BackgroundInnerRadiusVector[i] = 0.0;
419 BackgroundOuterRadiusVector[i] = 0.0;
420 SatellitePeakRadiusVector[i] = 0.0;
421 SatelliteBackgroundInnerRadiusVector[i] = 0.0;
422 SatelliteBackgroundOuterRadiusVector[i] = 0.0;
423 } else {
424 // Integrate peak properly
425 double inti;
426 double sigi;
427 std::vector<double> axes_radii;
428
429 // calculate adaptive background inner and outer radius
430 // compute adaptive background radius
431 double adaptiveBack_inner_radius = isSatellitePeak
432 ? adaptiveQBackgroundMultiplier * lenQpeak + satellite_back_inner_radius
433 : adaptiveQBackgroundMultiplier * lenQpeak + back_inner_radius;
434 double adaptiveBack_outer_radius = isSatellitePeak
435 ? adaptiveQBackgroundMultiplier * lenQpeak + satellite_back_outer_radius
436 : adaptiveQBackgroundMultiplier * lenQpeak + back_outer_radius;
437
438 // integrate the peak to get intensity and error
440 if (isSatellitePeak) {
441 // Satellite peak
442 SatellitePeakRadiusVector[i] = adaptiveRadius;
443 SatelliteBackgroundInnerRadiusVector[i] = adaptiveBack_inner_radius;
444 SatelliteBackgroundOuterRadiusVector[i] = adaptiveBack_outer_radius;
445
446 std::pair<double, double> backi;
448 if (!shareBackground || (satellitePeaks.size() > 0 &&
449 std::find(satellitePeaks.begin(), satellitePeaks.end(), i) != satellitePeaks.end())) {
450 // check if this satellite peak did NOT have a bragg peak, then we want to integrate it normally
451 shape =
452 integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, adaptiveRadius, adaptiveBack_inner_radius,
453 adaptiveBack_outer_radius, axes_radii, inti, sigi, backi);
454 } else {
455 // force satellite background radii in containers to use bragg peak background values
456 SatelliteBackgroundInnerRadiusVector[i] = adaptiveQBackgroundMultiplier * lenQpeak + back_inner_radius;
457 SatelliteBackgroundOuterRadiusVector[i] = adaptiveQBackgroundMultiplier * lenQpeak + back_outer_radius;
458
459 // if sharing background, integrate with background radii = peak radius so that background is 0 for now
460 shape = integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, adaptiveRadius, adaptiveRadius,
461 adaptiveRadius, axes_radii, inti, sigi, backi);
462 }
463
464 } else {
465 // Bragg peak
466 PeakRadiusVector[i] = adaptiveRadius;
467 BackgroundInnerRadiusVector[i] = adaptiveBack_inner_radius;
468 BackgroundOuterRadiusVector[i] = adaptiveBack_outer_radius;
469
470 std::pair<double, double> backi;
471 integrator.setRadius(m_braggPeakRadius);
472 shape =
473 integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, adaptiveRadius, adaptiveBack_inner_radius,
474 adaptiveBack_outer_radius, axes_radii, inti, sigi, backi);
475 if (shareBackground) {
476 // cache this bragg peak's background so we can apply it to all its satellite peaks later
477 cachedBraggBackground[i] = backi;
478 }
479 }
480
481 peaks[i].setIntensity(inti);
482 peaks[i].setSigmaIntensity(sigi);
483 peaks[i].setPeakShape(shape);
484 if (axes_radii.size() == 3) {
485 if (inti / sigi > cutoffIsigI || cutoffIsigI == EMPTY_DBL()) {
486 principalaxis1.emplace_back(axes_radii[0]);
487 principalaxis2.emplace_back(axes_radii[1]);
488 principalaxis3.emplace_back(axes_radii[2]);
489 }
490 }
491 }
492 }
493
494 // Remove background if backgrounds are shared
495 if (shareBackground) {
496 removeSharedBackground(satellitePeakMap, cachedBraggBackground);
497 }
498
499 if (principalaxis1.size() > 1) {
500 outputAxisProfiles(principalaxis1, principalaxis2, principalaxis3, cutoffIsigI, numSigmas, peaks, integrator);
501 }
502
503 // This flag is used by the PeaksWorkspace to evaluate whether it has been
504 // integrated.
505 peak_ws->mutableRun().addProperty("PeaksIntegrated", 1, true);
506 // These flags are specific to the algorithm.
507 peak_ws->mutableRun().addProperty("PeakRadius", PeakRadiusVector, true);
508 peak_ws->mutableRun().addProperty("BackgroundInnerRadius", BackgroundInnerRadiusVector, true);
509 peak_ws->mutableRun().addProperty("BackgroundOuterRadius", BackgroundOuterRadiusVector, true);
510 // These falgs are related to the satellite peaks and specific to the algorithm.
511 peak_ws->mutableRun().addProperty("SatellitePeakRadius", SatellitePeakRadiusVector, true);
512 peak_ws->mutableRun().addProperty("SatelliteBackgroundInnerRadius", SatelliteBackgroundInnerRadiusVector, true);
513 peak_ws->mutableRun().addProperty("SatelliteBackgroundOuterRadius", SatelliteBackgroundOuterRadiusVector, true);
514
515 setProperty("OutputWorkspace", peak_ws);
516}
517
519 m_targWSDescr.setMinMax(std::vector<double>(3, -2000.), std::vector<double>(3, 2000.));
522
523 // generate the detectors table
524 Mantid::API::Algorithm_sptr childAlg = createChildAlgorithm("PreprocessDetectorsToMD", 0.,
525 .5); // HACK. soft dependency on non-dependent package.
526 childAlg->setProperty("InputWorkspace", wksp);
527 childAlg->executeAsChildAlg();
528
529 DataObjects::TableWorkspace_sptr table = childAlg->getProperty("OutputWorkspace");
530 if (!table)
531 throw(std::runtime_error("Can not retrieve results of \"PreprocessDetectorsToMD\""));
532 else
534}
535
537 for (size_t i = 0; i < detectorInfo.size(); ++i) {
538 if (detectorInfo.isMonitor(i))
539 continue; // skip monitor
540 if (!detectorInfo.isMasked(i))
541 continue; // edge is masked so don't check if not masked
542 const auto &det = detectorInfo.detector(i);
543 double tt1 = det.getTwoTheta(V3D(0, 0, 0), V3D(0, 0, 1)); // two theta
544 double ph1 = det.getPhi(); // phi
545 V3D E1 = V3D(-std::sin(tt1) * std::cos(ph1), -std::sin(tt1) * std::sin(ph1),
546 1. - std::cos(tt1)); // end of trajectory
547 E1 = E1 * (1. / E1.norm()); // normalize
548 E1Vec.emplace_back(E1);
549 }
550}
551
555void IntegrateEllipsoidsV2::outputProfileWS(const std::vector<double> &principalaxis1,
556 const std::vector<double> &principalaxis2,
557 const std::vector<double> &principalaxis3, const std::string &wsname) {
558
559 constexpr size_t histogramNumber = 3;
560 Workspace_sptr wsProfile =
561 WorkspaceFactory::Instance().create("Workspace2D", histogramNumber, principalaxis1.size(), principalaxis1.size());
562 Workspace2D_sptr wsProfile2D = std::dynamic_pointer_cast<Workspace2D>(wsProfile);
563 AnalysisDataService::Instance().addOrReplace(wsname, wsProfile2D);
564
565 // set output workspace
566 Points points(principalaxis1.size(), LinearGenerator(0, 1));
567 wsProfile2D->setHistogram(0, points, Counts(std::move(principalaxis1)));
568 wsProfile2D->setHistogram(1, points, Counts(std::move(principalaxis2)));
569 wsProfile2D->setHistogram(2, points, Counts(std::move(principalaxis3)));
570}
571
579void IntegrateEllipsoidsV2::pairBraggSatellitePeaks(const size_t &n_peaks, std::vector<DataObjects::Peak> &peaks,
580 std::map<size_t, std::vector<Peak *>> &satellitePeakMap,
581 std::vector<size_t> &satellitePeaks) {
582
583 std::vector<size_t> braggPeaks;
584
585 for (size_t i = 0; i < n_peaks; i++) {
586 // check if peak is satellite peak
587 const bool isSatellitePeak = (peaks[i].getIntMNP().norm2() > 0);
588 // grab QLabFrame
589 const V3D peak_q = peaks[i].getQLabFrame();
590 // check if peak is origin (skip if true)
591 const bool isOrigin = isSatellitePeak ? IntegrateQLabEvents::isOrigin(peak_q, m_satellitePeakRadius)
593 if (isOrigin) {
594 continue;
595 }
596
597 if (isSatellitePeak) {
598 satellitePeaks.emplace_back(i);
599 } else {
600 braggPeaks.emplace_back(i);
601 }
602 }
603
604 // Generate mapping of all satellite peaks for each bragg peak
605 for (auto it = braggPeaks.begin(); it != braggPeaks.end(); it++) {
606 const auto braggHKL = peaks[*it].getIntHKL();
607
608 // loop over all satellite peaks to determine if it belongs to this bragg
609 for (auto satIt = satellitePeaks.begin(); satIt != satellitePeaks.end();) {
610 const auto satHKL = peaks[*satIt].getIntHKL();
611 if (satHKL == braggHKL) {
612 // this satellite peak shares the HKL vector, so it is a satellite peak of this bragg peak
613 satellitePeakMap[*it].emplace_back(&peaks[*satIt]);
614
615 // remove this sat peak from the list, since it can be associated with only one bragg peak
616 satIt = satellitePeaks.erase(satIt);
617 continue;
618 }
619 satIt++;
620 }
621 }
622
623 // Any leftover satellite peaks in this list means these did not have a bragg peak
624 if (satellitePeaks.size() > 0) {
625 g_log.debug() << "Unable to find Bragg peaks for " << satellitePeaks.size()
626 << " satellite peaks.. integrating these using the satellite background radii options.\n";
627 }
628}
629
633void IntegrateEllipsoidsV2::removeSharedBackground(std::map<size_t, std::vector<Peak *>> &satellitePeakMap,
634 std::map<size_t, std::pair<double, double>> &cachedBraggBackground) {
635
636 // loop over all bragg peaks and apply the cached background to their satellite peaks
637 for (auto it = satellitePeakMap.begin(); it != satellitePeakMap.end(); it++) {
638 const double bkgd_value{cachedBraggBackground[it->first].first};
639 const double bkgd_sigma{cachedBraggBackground[it->first].second};
640 for (auto satPeak = it->second.begin(); satPeak != it->second.end(); satPeak++) {
641 // subtract the cached background from the intensity
642 // (*satPeak)->setIntensity((*satPeak)->getIntensity() - cachedBraggBackground[it->first].first);
643 (*satPeak)->setIntensity((*satPeak)->getIntensity() - bkgd_value);
644
645 // update the sigma intensity based on the new background
646 double sigInt = (*satPeak)->getSigmaIntensity();
647 (*satPeak)->setSigmaIntensity(sqrt(sigInt * sigInt + bkgd_sigma));
648 }
649 }
650}
651
656void IntegrateEllipsoidsV2::outputAxisProfiles(std::vector<double> &principalaxis1, std::vector<double> &principalaxis2,
657 std::vector<double> &principalaxis3, const double &cutoffIsigI,
658 const int &numSigmas, std::vector<Peak> &peaks,
659 IntegrateQLabEvents &integrator) {
660
661 // Export principle axis profile to Fixed workspace EllipsoidAxes
662 outputProfileWS(principalaxis1, principalaxis2, principalaxis3, "EllipsoidAxes");
663
664 // Output message
665 Statistics stats1 = getStatistics(principalaxis1);
666 g_log.notice() << "principalaxis1: "
667 << " mean " << stats1.mean << " standard_deviation " << stats1.standard_deviation << " minimum "
668 << stats1.minimum << " maximum " << stats1.maximum << " median " << stats1.median << "\n";
669 Statistics stats2 = getStatistics(principalaxis2);
670 g_log.notice() << "principalaxis2: "
671 << " mean " << stats2.mean << " standard_deviation " << stats2.standard_deviation << " minimum "
672 << stats2.minimum << " maximum " << stats2.maximum << " median " << stats2.median << "\n";
673 Statistics stats3 = getStatistics(principalaxis3);
674 g_log.notice() << "principalaxis3: "
675 << " mean " << stats3.mean << " standard_deviation " << stats3.standard_deviation << " minimum "
676 << stats3.minimum << " maximum " << stats3.maximum << " median " << stats3.median << "\n";
677
678 // Some special case to amend ... ...
679 // Re-integrate peaks
680 if (cutoffIsigI != EMPTY_DBL()) {
681 double meanMax = std::max(std::max(stats1.mean, stats2.mean), stats3.mean);
682 double stdMax = std::max(std::max(stats1.standard_deviation, stats2.standard_deviation), stats3.standard_deviation);
683 integratePeaksCutoffISigI(meanMax, stdMax, principalaxis1, principalaxis2, principalaxis3, numSigmas, peaks,
684 integrator);
685
686 if (principalaxis1.size() > 1) {
687 outputProfileWS(principalaxis1, principalaxis2, principalaxis3, "EllipsoidAxes_2ndPass");
688 }
689 }
690}
691
696void IntegrateEllipsoidsV2::integratePeaksCutoffISigI(const double &meanMax, const double &stdMax,
697 std::vector<double> &principalaxis1,
698 std::vector<double> &principalaxis2,
699 std::vector<double> &principalaxis3, const int &numSigmas,
700 std::vector<Peak> &peaks, IntegrateQLabEvents &integrator) {
701 // reset all principle axes
702 principalaxis1.clear();
703 principalaxis2.clear();
704 principalaxis3.clear();
705
706 bool specify_size = true;
707 // double meanMax = std::max(std::max(stats1.mean, stats2.mean), stats3.mean);
708 // double stdMax = std::max(std::max(stats1.standard_deviation, stats2.standard_deviation),
709 // stats3.standard_deviation);
710 double peak_radius = meanMax + numSigmas * stdMax;
711 double back_inner_radius = peak_radius;
712 double back_outer_radius = peak_radius * 1.25992105; // A factor of 2 ^ (1/3)
713 // will make the background shell volume equal to the peak region volume.
714 for (size_t i = 0; i < peaks.size(); i++) {
715 // check if peak is satellite peak
716 const bool isSatellitePeak = (peaks[i].getIntMNP().norm2() > 0);
717 //
718 const V3D peak_q = peaks[i].getQLabFrame();
719 std::vector<double> axes_radii;
720
721 double inti{0.}, sigi{0.};
722 std::pair<double, double> backi;
723 if (isSatellitePeak) {
725 integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, peak_radius, back_inner_radius, back_outer_radius,
726 axes_radii, inti, sigi, backi);
727 } else {
728 integrator.setRadius(m_braggPeakRadius);
729 integrator.ellipseIntegrateEvents(E1Vec, peak_q, specify_size, peak_radius, back_inner_radius, back_outer_radius,
730 axes_radii, inti, sigi, backi);
731 }
732
733 peaks[i].setIntensity(inti);
734 peaks[i].setSigmaIntensity(sigi);
735 if (axes_radii.size() == 3) {
736 principalaxis1.emplace_back(axes_radii[0]);
737 principalaxis2.emplace_back(axes_radii[1]);
738 principalaxis3.emplace_back(axes_radii[2]);
739 }
740 }
741}
742
744 const std::string &property, const std::string &values) {
745 auto alg = createChildAlgorithm("MaskBTP");
746 alg->setProperty<Workspace_sptr>("Workspace", peakWS);
747 alg->setProperty(property, values);
748 if (!alg->execute())
749 throw std::runtime_error("MaskDetectors Child Algorithm has not executed successfully");
750}
751} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#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...
Definition: MultiThreaded.h:94
#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.
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
Kernel::Property * getPointerToProperty(const std::string &name) const override
Get a property by name.
Definition: Algorithm.cpp:2033
void removeProperty(const std::string &name, const bool delproperty=true) override
Removes the property from management.
Definition: Algorithm.cpp:2109
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
Definition: Algorithm.cpp:2076
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.
Definition: Algorithm.cpp:842
Kernel::Logger & g_log
Definition: Algorithm.h:451
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A class for holding :
Definition: EventList.h:56
void compressEvents(double tolerance, EventList *destination)
Compress the event list by grouping events with the same TOF (within a given tolerance).
Definition: EventList.cpp:1676
std::vector< WeightedEventNoTime > & getWeightedEventsNoTime()
Return the list of WeightedEvent contained.
Definition: EventList.cpp:823
void switchTo(Mantid::API::EventType newType) override
Switch the EventList to use the given EventType (TOF, WEIGHTED, or WEIGHTED_NOTIME)
Definition: EventList.cpp:649
bool empty() const
Much like stl containers, returns true if there is nothing in the event list.
Definition: EventList.cpp:1158
Geometry::DetectorInfo is an intermediate step towards a DetectorInfo that is part of Instrument-2....
Definition: DetectorInfo.h:49
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.
BoundedValidator is a validator that requires the values to be between upper or lower bounds,...
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
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 error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
virtual bool isDefault() const =0
Overriden function that returns if property has the same value that it was initialised with,...
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Class for 3D vectors.
Definition: V3D.h:34
double norm() const noexcept
Definition: V3D.h:263
void outputProfileWS(const std::vector< double > &principalaxis1, const std::vector< double > &principalaxis2, const std::vector< double > &principalaxis3, const std::string &wsname)
Write Axis profile to a MatrixWorkspace (Workspace2D)
void integratePeaksCutoffISigI(const double &meanMax, const double &stdMax, std::vector< double > &principalaxis1, std::vector< double > &principalaxis2, std::vector< double > &principalaxis3, const int &numSigmas, std::vector< DataObjects::Peak > &peaks, IntegrateQLabEvents &integrator_satellite)
Integrate peaks again with cutoff value of I/Sig(I)
void outputAxisProfiles(std::vector< double > &principalaxis1, std::vector< double > &principalaxis2, std::vector< double > &principalaxis3, const double &cutoffIsigI, const int &numSigmas, std::vector< DataObjects::Peak > &peaks, IntegrateQLabEvents &integrator)
Write the profiles of each principle axis to the output workspace with fixed name.
void removeSharedBackground(std::map< size_t, std::vector< DataObjects::Peak * > > &satellitePeakMap, std::map< size_t, std::pair< double, double > > &cachedBraggBackground)
Remove shared background from each satellite peak.
double m_satellitePeakRadius
peak radius for satellite peaks
std::vector< Kernel::V3D > E1Vec
save for all detector pixels
std::map< std::string, std::string > validateInputs() override
Private validator for inputs.
double m_braggPeakRadius
peak radius for Bragg peaks
void exec() override
Execute the algorithm.
void runMaskDetectors(const Mantid::DataObjects::PeaksWorkspace_sptr &peakWS, const std::string &property, const std::string &values)
void pairBraggSatellitePeaks(const size_t &n_peaks, std::vector< DataObjects::Peak > &peaks, std::map< size_t, std::vector< DataObjects::Peak * > > &satellitePeakMap, std::vector< size_t > &satellitePeaks)
Pair all Bragg peaks with their related satellite peaks.
void initTargetWSDescr(API::MatrixWorkspace_sptr &wksp)
Initialize the output information for the MD conversion framework.
void init() override
Initialize the algorithm's properties.
void qListFromEventWS(IntegrateQLabEvents &integrator, API::Progress &prog, DataObjects::EventWorkspace_sptr &wksp)
create a list of SlimEvent objects from an events workspace
void calculateE1(const Geometry::DetectorInfo &detectorInfo)
Calculate if this Q is on a detector.
void qListFromHistoWS(IntegrateQLabEvents &integrator, API::Progress &prog, DataObjects::Workspace2D_sptr &wksp)
create a list of SlimEvent objects from a histogram workspace
This is a low-level class to construct a map with lists of events near each peak Q-vector in the lab ...
static bool isOrigin(const V3D &q, const double &cellSize)
Determine if an input Q-vector lies in the cell associated to the origin.
void setRadius(const double &radius)
Set peak integration radius.
void populateCellsWithPeaks()
Assign events to each of the cells occupied by events.
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, std::pair< double, double > &backi)
Integrate the events around the specified peak QLab vector.
void addEvents(SlimEvents const &event_qs)
distribute the events among the cells of the partitioned QLab space.
Class responsible for conversion of input workspace data into proper number of output dimensions for ...
Definition: MDTransfQ3D.h:31
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...
Definition: MDTransfQ3D.cpp:47
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.
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
Definition: Workspace_fwd.h:20
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
Definition: Algorithm.h:61
@ WEIGHTED_NOTIME
Definition: IEventList.h:18
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::shared_ptr< const PeakShape > PeakShape_const_sptr
Definition: PeakShape.h:43
Statistics getStatistics(const std::vector< TYPE > &data, const unsigned int flags=StatOptions::AllStats)
Return a statistics object for the given data set.
Definition: Statistics.cpp:167
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.
Definition: MultiThreaded.h:22
std::vector< SlimEvent > SlimEvents
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.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
Simple struct to store statistics.
Definition: Statistics.h:25
double mean
Mean value.
Definition: Statistics.h:31
double median
Median value.
Definition: Statistics.h:33
double minimum
Minimum value.
Definition: Statistics.h:27
double maximum
Maximum value.
Definition: Statistics.h:29
double standard_deviation
standard_deviation of the values
Definition: Statistics.h:35