Mantid
Loading...
Searching...
No Matches
DgsDiagnose.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 +
12
13#include <boost/lexical_cast.hpp>
14#include <boost/pointer_cast.hpp>
15
16using namespace Mantid::Kernel;
17using namespace Mantid::API;
18using namespace Mantid::DataObjects;
19using namespace WorkflowAlgorithmHelpers;
20
22
23// Register the algorithm into the AlgorithmFactory
24DECLARE_ALGORITHM(DgsDiagnose)
25
26//----------------------------------------------------------------------------------------------
28const std::string DgsDiagnose::name() const { return "DgsDiagnose"; }
29
31int DgsDiagnose::version() const { return 1; }
32
34const std::string DgsDiagnose::category() const { return "Workflow\\Inelastic\\UsesPropertyManager"; }
35
36//----------------------------------------------------------------------------------------------
37
38//----------------------------------------------------------------------------------------------
42 this->declareProperty(std::make_unique<WorkspaceProperty<>>("DetVanWorkspace", "", Direction::Input),
43 "The detector vanadium workspace.");
44 this->declareProperty(
45 std::make_unique<WorkspaceProperty<>>("DetVanMonitorWorkspace", "", Direction::Input, PropertyMode::Optional),
46 "A monitor workspace associated with the detector vanadium workspace.");
47 this->declareProperty(
48 std::make_unique<WorkspaceProperty<>>("DetVanCompWorkspace", "", Direction::Input, PropertyMode::Optional),
49 "A detector vanadium workspace to compare against the primary one.");
50 this->declareProperty(
51 std::make_unique<WorkspaceProperty<>>("DetVanCompMonitorWorkspace", "", Direction::Input, PropertyMode::Optional),
52 "A monitor workspace associated with the comparison "
53 "detector vanadium workspace.");
54 this->declareProperty(
55 std::make_unique<WorkspaceProperty<>>("SampleWorkspace", "", Direction::Input, PropertyMode::Optional),
56 "A sample workspace to run some diagnostics on.");
57 this->declareProperty(
58 std::make_unique<WorkspaceProperty<>>("SampleMonitorWorkspace", "", Direction::Input, PropertyMode::Optional),
59 "A monitor workspace associated with the sample workspace.");
60 this->declareProperty(
61 std::make_unique<WorkspaceProperty<>>("HardMaskWorkspace", "", Direction::Input, PropertyMode::Optional),
62 "A hard mask workspace to apply.");
63 this->declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
64 "This is the resulting mask workspace.");
65 this->declareProperty("ReductionProperties", "__dgs_reduction_properties", Direction::Input);
66}
67
68//----------------------------------------------------------------------------------------------
72 g_log.notice() << "Starting DgsDiagnose\n";
73 // Get the reduction property manager
74 const std::string reductionManagerName = this->getProperty("ReductionProperties");
75 std::shared_ptr<PropertyManager> reductionManager;
76 if (PropertyManagerDataService::Instance().doesExist(reductionManagerName)) {
77 reductionManager = PropertyManagerDataService::Instance().retrieve(reductionManagerName);
78 } else {
79 throw std::runtime_error("DgsDiagnose cannot run without a reduction PropertyManager.");
80 }
81
82 // Gather all the necessary properties
83 MatrixWorkspace_sptr detVanWS = this->getProperty("DetVanWorkspace");
84 MatrixWorkspace_sptr detVanMonWS = this->getProperty("DetVanMonitorWorkspace");
85 MatrixWorkspace_sptr detVanCompWS = this->getProperty("DetVanCompWorkspace");
86 MatrixWorkspace_sptr detVanCompMonWS = this->getProperty("DetVanCompMonitorWorkspace");
87 MatrixWorkspace_sptr hardMaskWS = this->getProperty("HardMaskWorkspace");
88 MatrixWorkspace_sptr sampleWS;
89 MatrixWorkspace_sptr sampleMonWS;
90
91 // Boolean properties
92 const bool checkBkg = getBoolPropOrParam("BackgroundCheck", reductionManager, "check_background", detVanWS);
93 const bool rejectZeroBkg = getBoolPropOrParam("RejectZeroBackground", reductionManager, "diag_samp_zero", detVanWS);
94 const bool createPsdBleed = getBoolPropOrParam("PsdBleed", reductionManager, "diag_bleed_test", detVanWS);
95 const bool vanSA =
96 getBoolPropOrParam("MedianTestCorrectForSolidAngle", reductionManager, "diag_correct_solid_angle", detVanWS);
97
98 // Numeric properties
99 const double huge = getDblPropOrParam("HighCounts", reductionManager, "diag_huge", detVanWS);
100 const double tiny = getDblPropOrParam("LowCounts", reductionManager, "diag_tiny", detVanWS);
101 const double vanOutHi = getDblPropOrParam("HighOutlier", reductionManager, "diag_van_out_hi", detVanWS);
102 const double vanOutLo = getDblPropOrParam("LowOutlier", reductionManager, "diag_van_out_lo", detVanWS);
103 const double vanHi = getDblPropOrParam("MedianTestHigh", reductionManager, "diag_van_hi", detVanWS);
104 const double vanLo = getDblPropOrParam("MedianTestLow", reductionManager, "diag_van_lo", detVanWS);
105 const double vanLevelsUp = getDblPropOrParam("MedianTestLevelsUp", reductionManager, "diag_van_levels", detVanWS, 0);
106 const double vanSigma = getDblPropOrParam("ErrorBarCriterion", reductionManager, "diag_van_sig", detVanWS);
107 const double variation = getDblPropOrParam("DetVanRatioVariation", reductionManager, "diag_variation", detVanWS);
108 const double samHi = getDblPropOrParam("SamBkgMedianTestHigh", reductionManager, "diag_samp_hi", detVanWS);
109 const double samLo = getDblPropOrParam("SamBkgMedianTestLow", reductionManager, "diag_samp_lo", detVanWS);
110 const double samSigma = getDblPropOrParam("SamBkgErrorBarCriterion", reductionManager, "diag_samp_sig", detVanWS);
111 double bleedRate = getDblPropOrParam("MaxFramerate", reductionManager, "diag_bleed_maxrate", detVanWS);
112 const double bleedPixels = getDblPropOrParam("IgnoredPixels", reductionManager, "diag_bleed_pixels", detVanWS, 80.0);
113
114 // Make some internal names for workspaces
115 const std::string dvInternal = "__det_van";
116 const std::string dvCompInternal = "__det_van_comp";
117 const std::string sampleInternal = "__sample";
118 const std::string bkgInternal = "__background_int";
119 const std::string countsInternal = "__total_counts";
120
121 // If we are running this standalone, the IncidentEnergyGuess property in
122 // the reduction property manager does not exist. If that is true, then we
123 // don't have to clone workspaces.
124 bool isStandAlone = !reductionManager->existsProperty("IncidentEnergyGuess");
125
126 // Process the detector vanadium
127 auto detVan = createChildAlgorithm("DgsProcessDetectorVanadium");
128 detVan->setProperty("InputWorkspace", detVanWS);
129 detVan->setProperty("OutputWorkspace", dvInternal);
130 detVan->setProperty("InputMonitorWorkspace", detVanMonWS);
131 detVan->setProperty("ReductionProperties", reductionManagerName);
132 detVan->executeAsChildAlg();
133 MatrixWorkspace_sptr dvWS = detVan->getProperty("OutputWorkspace");
134
135 // Process the comparison detector vanadium workspace if present
136 MatrixWorkspace_sptr dvCompWS;
137 if (detVanCompWS) {
138 detVan->setProperty("InputWorkspace", detVanCompWS);
139 detVan->setProperty("OutputWorkspace", dvCompInternal);
140 detVan->setProperty("InputMonitorWorkspace", detVanCompMonWS);
141 detVan->executeAsChildAlg();
142 dvCompWS = detVan->getProperty("OutputWorkspace");
143 detVanCompWS.reset();
144 } else {
145 dvCompWS = std::shared_ptr<MatrixWorkspace>();
146 }
147
148 // Process the sample data if any of the sample checks are requested.
149 if (checkBkg || rejectZeroBkg || createPsdBleed) {
150 sampleWS = this->getProperty("SampleWorkspace");
151 sampleMonWS = this->getProperty("SampleMonitorWorkspace");
152
154 if (!isStandAlone) {
155 auto cloneWs = createChildAlgorithm("CloneWorkspace");
156 cloneWs->setProperty("InputWorkspace", sampleWS);
157 cloneWs->setProperty("OutputWorkspace", sampleInternal);
158 cloneWs->executeAsChildAlg();
159 tmp = cloneWs->getProperty("OutputWorkspace");
160 sampleWS = std::static_pointer_cast<MatrixWorkspace>(tmp);
161 }
162
163 auto norm = createChildAlgorithm("DgsPreprocessData");
164 norm->setProperty("InputWorkspace", sampleWS);
165 norm->setProperty("OutputWorkspace", sampleWS);
166 norm->setProperty("InputMonitorWorkspace", sampleMonWS);
167 norm->setProperty("ReductionProperties", reductionManagerName);
168 norm->executeAsChildAlg();
169 sampleWS = norm->getProperty("OutputWorkspace");
170 }
171
172 // Create the total counts workspace if necessary
173 MatrixWorkspace_sptr totalCountsWS;
174 if (rejectZeroBkg) {
175 auto integrate = createChildAlgorithm("Integration");
176 integrate->setProperty("InputWorkspace", sampleWS);
177 integrate->setProperty("OutputWorkspace", countsInternal);
178 integrate->setProperty("IncludePartialBins", true);
179 integrate->executeAsChildAlg();
180 totalCountsWS = integrate->getProperty("OutputWorkspace");
181 } else {
182 totalCountsWS = std::shared_ptr<MatrixWorkspace>();
183 }
184
185 // Create the background workspace if necessary
186 MatrixWorkspace_sptr backgroundIntWS;
187 if (checkBkg) {
188 double rangeStart = getDblPropOrParam("BackgroundTofStart", reductionManager, "bkgd-range-min", sampleWS);
189
190 double rangeEnd = getDblPropOrParam("BackgroundTofEnd", reductionManager, "bkgd-range-max", sampleWS);
191
192 auto integrate = createChildAlgorithm("Integration");
193 integrate->setProperty("InputWorkspace", sampleWS);
194 integrate->setProperty("OutputWorkspace", bkgInternal);
195 integrate->setProperty("RangeLower", rangeStart);
196 integrate->setProperty("RangeUpper", rangeEnd);
197 integrate->setProperty("IncludePartialBins", true);
198 integrate->executeAsChildAlg();
199 backgroundIntWS = integrate->getProperty("OutputWorkspace");
200
201 // Need to match the units between background and detector vanadium
202 const std::string detVanIntRangeUnits = reductionManager->getProperty("DetVanIntRangeUnits");
203 auto cvu = createChildAlgorithm("ConvertUnits");
204 cvu->setProperty("InputWorkspace", backgroundIntWS);
205 cvu->setProperty("OutputWorkspace", backgroundIntWS);
206 cvu->setProperty("Target", detVanIntRangeUnits);
207 cvu->executeAsChildAlg();
208 backgroundIntWS = cvu->getProperty("OutputWorkspace");
209
210 // Normalise the background integral workspace
211 if (dvCompWS) {
212 MatrixWorkspace_sptr hmean = 2.0 * dvWS * dvCompWS;
213 hmean /= (dvWS + dvCompWS);
214 backgroundIntWS /= hmean;
215 } else {
216 backgroundIntWS /= dvWS;
217 }
218 } else {
219 backgroundIntWS = std::shared_ptr<MatrixWorkspace>();
220 }
221
222 // Handle case where one of the other tests (checkBkg or rejectZeroBkg)
223 // are requested, but not createPsdBleed.
224 if (!createPsdBleed) {
225 sampleWS = std::shared_ptr<MatrixWorkspace>();
226 }
227
228 auto diag = createChildAlgorithm("DetectorDiagnostic");
229 diag->setProperty("InputWorkspace", dvWS);
230 diag->setProperty("DetVanCompare", dvCompWS);
231 diag->setProperty("SampleWorkspace", sampleWS);
232 diag->setProperty("SampleTotalCountsWorkspace", totalCountsWS);
233 diag->setProperty("SampleBackgroundWorkspace", backgroundIntWS);
234 diag->setProperty("HardMaskWorkspace", hardMaskWS);
235 diag->setProperty("LowThreshold", tiny);
236 diag->setProperty("HighThreshold", huge);
237 diag->setProperty("LowOutlier", vanOutLo);
238 diag->setProperty("HighOutlier", vanOutHi);
239 diag->setProperty("LowThresholdFraction", vanLo);
240 diag->setProperty("HighThresholdFraction", vanHi);
241 diag->setProperty("LevelsUp", static_cast<int>(vanLevelsUp));
242 diag->setProperty("CorrectForSolidAngle", vanSA);
243 diag->setProperty("SignificanceTest", vanSigma);
244 diag->setProperty("DetVanRatioVariation", variation);
245 diag->setProperty("SampleBkgLowAcceptanceFactor", samLo);
246 diag->setProperty("SampleBkgHighAcceptanceFactor", samHi);
247 diag->setProperty("SampleBkgSignificanceTest", samSigma);
248 diag->setProperty("MaxTubeFramerate", bleedRate);
249 diag->setProperty("NIgnoredCentralPixels", static_cast<int>(bleedPixels));
250
252 std::vector<std::string> diag_spectra = dvWS->getInstrument()->getStringParameter("diag_spectra");
253 if (diag_spectra.empty() || "None" == diag_spectra[0]) {
254 diag->execute();
255 maskWS = diag->getProperty("OutputWorkspace");
256 } else {
258 tokenizer tokens(diag_spectra[0], "(,);", Mantid::Kernel::StringTokenizer::TOK_IGNORE_EMPTY);
259 for (auto tok_iter = tokens.begin(); tok_iter != tokens.end();) {
260 auto startIndex = boost::lexical_cast<int>(*tok_iter);
261 startIndex -= 1;
262 ++tok_iter;
263 auto endIndex = boost::lexical_cast<int>(*tok_iter);
264 endIndex -= 1;
265 g_log.information() << "Pixel range: (" << startIndex << ", " << endIndex << ")\n";
266 diag->setProperty("StartWorkspaceIndex", startIndex);
267 diag->setProperty("EndWorkspaceIndex", endIndex);
268 diag->execute();
269 if (maskWS) {
270 MatrixWorkspace_sptr tmp = diag->getProperty("OutputWorkspace");
271 auto comb = createChildAlgorithm("BinaryOperateMasks");
272 comb->setProperty("InputWorkspace1", maskWS);
273 comb->setProperty("InputWorkspace2", tmp);
274 comb->setProperty("OutputWorkspace", maskWS);
275 comb->setProperty("OperationType", "OR");
276 comb->execute();
277 } else {
278 maskWS = diag->getProperty("OutputWorkspace");
279 }
280 ++tok_iter;
281 }
282 }
283
284 // Cleanup
285 dvWS.reset();
286 dvCompWS.reset();
287 sampleWS.reset();
288 totalCountsWS.reset();
289 backgroundIntWS.reset();
290
291 // If mask file name is set save out the diagnostic mask.
292 if (reductionManager->existsProperty("OutputMaskFile")) {
293 std::string maskFilename = reductionManager->getPropertyValue("OutputMaskFile");
294 if (!maskFilename.empty()) {
295 auto saveNxs = createChildAlgorithm("SaveMask");
296 saveNxs->setProperty("InputWorkspace", maskWS);
297 saveNxs->setProperty("OutputFile", maskFilename);
298 saveNxs->execute();
299 }
300 }
301
302 MaskWorkspace_sptr m = std::dynamic_pointer_cast<MaskWorkspace>(maskWS);
303 g_log.information() << "Number of masked pixels = " << m->getNumberMasked() << '\n';
304
305 this->setProperty("OutputWorkspace", maskWS);
306}
307
308} // namespace Mantid::WorkflowAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_vector * tmp
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
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
A property class for workspaces.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
Iterator begin()
Iterator referring to first element in the container.
@ TOK_IGNORE_EMPTY
ignore empty tokens
Iterator end()
Iterator referring to the past-the-end element in the container.
DgsDiagnose : This algorithm constructs all of the necessary workspaces for performing detector diagn...
Definition: DgsDiagnose.h:19
const std::string category() const override
Algorithm's category for identification.
Definition: DgsDiagnose.cpp:34
void init() override
Initialize the algorithm's properties.
Definition: DgsDiagnose.cpp:41
int version() const override
Algorithm's version for identification.
Definition: DgsDiagnose.cpp:31
void exec() override
Execute the algorithm.
Definition: DgsDiagnose.cpp:71
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< MaskWorkspace > MaskWorkspace_sptr
shared pointer to the MaskWorkspace class
Definition: MaskWorkspace.h:64
bool getBoolPropOrParam(const std::string &pmProp, Mantid::Kernel::PropertyManager_sptr &pm, const std::string &instParam, Mantid::API::MatrixWorkspace_sptr &ws, const bool overrideValue=false)
Function to get boolean property or instrument parameter value.
double getDblPropOrParam(const std::string &pmProp, Mantid::Kernel::PropertyManager_sptr &pm, const std::string &instParam, Mantid::API::MatrixWorkspace_sptr &ws, const double overrideValue=Mantid::EMPTY_DBL())
Function to get double property or instrument parameter value.
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54