Mantid
Loading...
Searching...
No Matches
PreprocessDetectorsToMD.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#include "MantidAPI/Run.h"
17
18using namespace Mantid;
19using namespace Mantid::API;
20using namespace Mantid::DataObjects;
21
22namespace Mantid::MDAlgorithms {
23// Register the algorithm into the AlgorithmFactory
24DECLARE_ALGORITHM(PreprocessDetectorsToMD)
25
26PreprocessDetectorsToMD::PreprocessDetectorsToMD() : m_getEFixed(false), m_getIsMasked(false) {}
27
30 auto ws_valid = std::make_shared<Kernel::CompositeValidator>();
31 // workspace needs instrument
32 ws_valid->add<API::InstrumentValidator>();
33 // the validator which checks if the workspace has axis and any units
34 // ws_valid->add<API::WorkspaceUnitValidator>("");
35
36 declareProperty(std::make_unique<WorkspaceProperty<API::MatrixWorkspace>>("InputWorkspace", "",
37 Kernel::Direction::Input, ws_valid),
38 "Name of an input Matrix Workspace with instrument.");
39
40 declareProperty(std::make_unique<WorkspaceProperty<TableWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
41 "Name of the output Table workspace with pre-processed "
42 "detectors data. If the workspace exists, it will be "
43 "replaced.");
44
45 declareProperty(std::make_unique<Kernel::PropertyWithValue<bool>>("GetMaskState", true, Kernel::Direction::Input),
46 "Returns masked state of the detectors. If this option is "
47 "false, the masked detectors are just dropped from the "
48 "resulting workspace and spectra-to detectors map has to be "
49 "used to analyse the spectra. This is temporary parameter "
50 "and logic necessary until Mantid masks signal by 0 rather "
51 "then NaN.");
52
53 declareProperty(std::make_unique<Kernel::PropertyWithValue<bool>>("UpdateMasksInfo", false, Kernel::Direction::Input),
54 "If target workspace already exists as the result of "
55 "previous deployment of this algorithm, the algorithm just "
56 "updated masks states column instead of calculating the "
57 "whole target workspace. The target workspace has to be "
58 "appropriate for the source workspace This is temporary "
59 "parameter and logic necessary until Mantid masks signal by "
60 "0 rather then NaN.");
61
63 "This option makes sense for Indirect instrument, where each "
64 "detector can have its own energy, defined by correspondent "
65 "crystal-analyzer position.\n"
66 "If this option is selected for other instrument types, the "
67 "value of eFixed is taken from workspace property "
68 "Ei"
69 " or "
70 "eFixed"
71 " if "
72 "Ei"
73 "\n"
74 "is missing and is set to NaN if no such properties are "
75 "defined on the input workspace.");
76
77 // declareProperty(new
78 // PropertyWithValue<bool>("FakeDetectors",false,Kernel::Direction::Input),
79 // "If selected, generates table workspace with fake detectors, all allocated
80 // in a monitor position.\n"
81 // "Number of detectors is equal to number of spectra and real detectors, if
82 // present are ignored ");
83}
84
85//----------------------------------------------------------------------------------------------
86/* Execute the algorithm. */
88 // -------- get Input workspace
89 MatrixWorkspace_const_sptr inputWS = getProperty("InputWorkspace");
90
91 // verify if we want to preprocess real workspace or just update the state of
92 // workspace masks
94 bool updateMasks(false);
95 if (this->getProperty("GetMaskState") && this->getProperty("UpdateMasksInfo")) {
96 std::string wsName = this->getPointerToProperty("OutputWorkspace")->value();
97 if (API::AnalysisDataService::Instance().doesExist(wsName)) {
98 targWS =
99 std::dynamic_pointer_cast<DataObjects::TableWorkspace>(API::AnalysisDataService::Instance().retrieve(wsName));
100 if (targWS) {
101 auto *pMasksArray = targWS->getColDataArray<int>("detMask");
102 if (pMasksArray)
103 updateMasks = true;
104 // was this workspace calculated without eFixed and now we need one?
105 if (this->getProperty("GetEFixed") && !targWS->getColDataArray<float>("eFixed"))
106 updateMasks = false;
107 }
108 }
109 }
110
111 if (updateMasks) // just update masks
112 this->updateMasksState(inputWS, targWS);
113 else // -------- build target workspace:
114 targWS = createTableWorkspace(inputWS);
115
116 if (this->isDetInfoLost(inputWS))
117 this->buildFakeDetectorsPositions(inputWS, targWS);
118 else // process real detectors positions
119 this->processDetectorsPositions(inputWS, targWS);
120
121 // set up target workspace
122 setProperty("OutputWorkspace", targWS);
123}
124
126std::shared_ptr<DataObjects::TableWorkspace>
128 const size_t nHist = inputWS->getNumberHistograms();
129
130 // set the target workspace
131 auto targWS = std::make_shared<TableWorkspace>(nHist);
132 // detectors positions
133 targWS->addColumn("V3D", "DetDirections");
134 // sample-detector distance;
135 targWS->addColumn("double", "L2");
136 // Diffraction angle
137 targWS->addColumn("double", "TwoTheta");
138 targWS->addColumn("double", "Azimuthal");
139 // the detector ID;
140 targWS->addColumn("int", "DetectorID");
141 // stores spectra index which corresponds to a valid detector index;
142 targWS->addColumn("size_t", "detIDMap");
143 // stores detector index which corresponds to the workspace index;
144 targWS->addColumn("size_t", "spec2detMap");
145
146 m_getIsMasked = this->getProperty("GetMaskState");
147 if (m_getIsMasked) // as bool is presented in vectors as a class, we are using
148 // int instead of bool
149 targWS->addColumn("int", "detMask");
150
151 // check if one wants to obtain detector's efixed"
152 m_getEFixed = this->getProperty("GetEFixed");
153 if (m_getEFixed)
154 targWS->addColumn("float", "eFixed");
155 targWS->addColumn("double", "DIFA");
156 targWS->addColumn("double", "DIFC");
157 targWS->addColumn("double", "TZERO");
158
159 // will see about that
160 // sin^2(Theta)
161 // std::vector<double> SinThetaSq;
162
163 double Efi = getEi(inputWS);
164 targWS->logs()->addProperty<double>("Ei", Efi, true);
165
166 return targWS;
167}
168
175 g_log.information() << "Preprocessing detector locations in a target reciprocal space\n";
176 //
177 Geometry::Instrument_const_sptr instrument = inputWS->getInstrument();
178 // this->pBaseInstr = instrument->baseInstrument();
179 //
180 Geometry::IComponent_const_sptr source = instrument->getSource();
181 Geometry::IComponent_const_sptr sample = instrument->getSample();
182 if ((!source) || (!sample)) {
183 g_log.error() << " Instrument is not fully defined. Can not identify "
184 "source or sample\n";
186 "Instrument not sufficiently defined: failed to get source and/or "
187 "sample");
188 }
189
190 // L1
191 try {
192 double L1 = source->getDistance(*sample);
193 targWS->logs()->addProperty<double>("L1", L1, true);
194 g_log.debug() << "Source-sample distance: " << L1 << '\n';
196 throw Kernel::Exception::InstrumentDefinitionError("Unable to calculate source-sample distance for workspace",
197 inputWS->getTitle());
198 }
199 // Instrument name
200 std::string InstrName = instrument->getName();
201 targWS->logs()->addProperty<std::string>("InstrumentName", InstrName,
202 true); // "The name which should unique identify current instrument");
203 targWS->logs()->addProperty<bool>("FakeDetectors", false, true);
204
205 // get access to the workspace memory
206 auto &sp2detMap = targWS->getColVector<size_t>("spec2detMap");
207 auto &detId = targWS->getColVector<int32_t>("DetectorID");
208 auto &detIDMap = targWS->getColVector<size_t>("detIDMap");
209 auto &L2 = targWS->getColVector<double>("L2");
210 auto &TwoTheta = targWS->getColVector<double>("TwoTheta");
211 auto &Azimuthal = targWS->getColVector<double>("Azimuthal");
212 auto &detDir = targWS->getColVector<Kernel::V3D>("DetDirections");
213 auto &DIFA = targWS->getColVector<double>("DIFA");
214 auto &DIFC = targWS->getColVector<double>("DIFC");
215 auto &TZERO = targWS->getColVector<double>("TZERO");
216
217 // Efixed; do we need one and does one exist?
218 auto Efi = targWS->getLogs()->getPropertyValueAsType<double>("Ei");
219 float *pEfixedArray(nullptr);
220 const Geometry::ParameterMap &pmap = inputWS->constInstrumentParameters();
221 if (m_getEFixed)
222 pEfixedArray = targWS->getColDataArray<float>("eFixed");
223
224 // check if one needs to generate masked detectors column.
225 int *pMasksArray(nullptr);
226 if (m_getIsMasked)
227 pMasksArray = targWS->getColDataArray<int>("detMask");
228
230 size_t div = 100;
231 size_t nHist = targWS->rowCount();
232 Mantid::API::Progress theProgress(this, 0.0, 1.0, nHist);
234 uint32_t liveDetectorsCount(0);
235 const auto &spectrumInfo = inputWS->spectrumInfo();
236 for (size_t i = 0; i < nHist; i++) {
237 sp2detMap[i] = std::numeric_limits<uint64_t>::quiet_NaN();
238 detId[i] = std::numeric_limits<int32_t>::quiet_NaN();
239 detIDMap[i] = std::numeric_limits<uint64_t>::quiet_NaN();
240 L2[i] = std::numeric_limits<double>::quiet_NaN();
241 DIFC[i] = std::numeric_limits<double>::quiet_NaN();
242 TwoTheta[i] = std::numeric_limits<double>::quiet_NaN();
243 Azimuthal[i] = std::numeric_limits<double>::quiet_NaN();
244 // detMask[i] = true;
245
246 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
247 continue;
248
249 // if masked detectors state is not used, masked detectors just ignored;
250 bool maskDetector = spectrumInfo.isMasked(i);
251 if (m_getIsMasked)
252 *(pMasksArray + liveDetectorsCount) = maskDetector ? 1 : 0;
253 else if (maskDetector)
254 continue;
255
256 const auto &spDet = spectrumInfo.detector(i);
257
258 // calculate the requested values;
259 sp2detMap[i] = liveDetectorsCount;
260 detId[liveDetectorsCount] = int32_t(spDet.getID());
261 detIDMap[liveDetectorsCount] = i;
262 L2[liveDetectorsCount] = spectrumInfo.l2(i);
263
264 double polar = spectrumInfo.twoTheta(i);
265 double azim = spDet.getPhi();
266 TwoTheta[liveDetectorsCount] = polar;
267 Azimuthal[liveDetectorsCount] = azim;
268
269 std::vector<int> warningDets;
270 auto diffConsts = spectrumInfo.diffractometerConstants(i, warningDets);
271 // map will create an entry with zero value if not present already
272 DIFA[liveDetectorsCount] = diffConsts[Kernel::UnitParams::difa];
273 DIFC[liveDetectorsCount] = diffConsts[Kernel::UnitParams::difc];
274 TZERO[liveDetectorsCount] = diffConsts[Kernel::UnitParams::tzero];
275
276 double sPhi = sin(polar);
277 double ez = cos(polar);
278 double ex = sPhi * cos(azim);
279 double ey = sPhi * sin(azim);
280
281 detDir[liveDetectorsCount].setX(ex);
282 detDir[liveDetectorsCount].setY(ey);
283 detDir[liveDetectorsCount].setZ(ez);
284
285 // double sinTheta=sin(0.5*polar);
286 // this->SinThetaSq[liveDetectorsCount] = sinTheta*sinTheta;
287
288 // specific code which should work and makes sense
289 // for indirect instrument but may be deployed on any code with Ei property
290 // defined;
291 if (pEfixedArray) {
292 try {
293 Geometry::Parameter_sptr par = pmap.getRecursive(&spDet, "eFixed");
294 if (par)
295 Efi = par->value<double>();
296 } catch (std::runtime_error &) {
297 }
298 // set efixed for each existing detector
299 *(pEfixedArray + liveDetectorsCount) = static_cast<float>(Efi);
300 }
301
302 liveDetectorsCount++;
303 if (i % div == 0)
304 theProgress.report(i, "Preprocessing detectors");
305 }
306 targWS->logs()->addProperty<uint32_t>("ActualDetectorsNum", liveDetectorsCount, true);
307
308 theProgress.report();
309 g_log.information() << "Finished preprocessing detector locations. Found: " << liveDetectorsCount
310 << " detectors out of: " << nHist << " histograms\n";
311}
318 auto *pMasksArray = targWS->getColDataArray<int>("detMask");
319 if (!pMasksArray)
320 throw std::invalid_argument("target workspace " + targWS->getName() +
321 " does not have defined masks column to update");
322
323 size_t nHist = targWS->rowCount();
324 const size_t nRows = inputWS->getNumberHistograms();
325 if (nHist != nRows)
326 throw std::invalid_argument(" source workspace " + inputWS->getName() + " and target workspace " +
327 targWS->getName() + " are inconsistent as have different numner of detectors");
328
329 uint32_t liveDetectorsCount(0);
330 const auto &spectrumInfo = inputWS->spectrumInfo();
331 for (size_t i = 0; i < nHist; i++) {
332 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
333 continue;
334
335 // if masked detectors state is not used, masked detectors just ignored;
336 bool maskDetector = spectrumInfo.isMasked(i);
337 *(pMasksArray + liveDetectorsCount) = maskDetector ? 1 : 0;
338
339 liveDetectorsCount++;
340 }
341}
342
347 UNUSED_ARG(inputWS);
348 // set sample-detector position equal to 1;
349 targWS->logs()->addProperty<double>("L1", 1., true);
350 //
351 targWS->logs()->addProperty<std::string>("InstrumentName", "FakeInstrument", true);
352 targWS->logs()->addProperty<bool>("FakeDetectors", true, true);
353
354 // get access to the workspace memory
355 auto &sp2detMap = targWS->getColVector<size_t>("spec2detMap");
356 auto &detId = targWS->getColVector<int32_t>("DetectorID");
357 auto &detIDMap = targWS->getColVector<size_t>("detIDMap");
358 auto &L2 = targWS->getColVector<double>("L2");
359 auto &TwoTheta = targWS->getColVector<double>("TwoTheta");
360 auto &Azimuthal = targWS->getColVector<double>("Azimuthal");
361 auto &detDir = targWS->getColVector<Kernel::V3D>("DetDirections");
362 // auto &detMask = targWS->getColVector<bool>("detMask");
363
365 size_t nHist = targWS->rowCount();
366 targWS->logs()->addProperty<uint32_t>("ActualDetectorsNum", uint32_t(nHist), true);
367
368 double polar(0);
369 // Loop over the spectra
370 for (size_t i = 0; i < nHist; i++) {
371 sp2detMap[i] = i;
372 detId[i] = static_cast<detid_t>(i);
373 detIDMap[i] = i;
374 L2[i] = 1;
375
376 TwoTheta[i] = polar;
377 Azimuthal[i] = 0;
378 // this->SinThetaSq[i]= 0;
379
380 double ez = 1.;
381 double ex = 0.;
382 double ey = 0.;
383
384 detDir[i].setX(ex);
385 detDir[i].setY(ey);
386 detDir[i].setZ(ez);
387 }
388 //
389}
390
394 auto pYAxis = dynamic_cast<API::NumericAxis *>(inWS2D->getAxis(1));
395 // if this is numeric axis, then the detector's information has been lost:
396 return pYAxis != nullptr;
397}
398
412 double Efi = std::numeric_limits<double>::quiet_NaN();
413
414 // is Ei on workspace properties? (it can be defined for some reason if
415 // detectors do not have one, and then it would exist as Ei)
416 bool EiFound(false);
417 try {
418 Efi = inputWS->run().getPropertyValueAsType<double>("Ei");
419 EiFound = true;
421 }
422 // try to get Efixed as property on a workspace, obtained for indirect
423 // instrument
424 bool eFixedFound(false);
425 if (!EiFound) {
426 try {
427 Efi = inputWS->run().getPropertyValueAsType<double>("eFixed");
428 eFixedFound = true;
430 }
431 }
432
433 if (!(EiFound || eFixedFound))
434 g_log.debug() << " Ei/eFixed requested but have not been found\n";
435
436 return Efi;
437}
438} // namespace Mantid::MDAlgorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
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
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
A validator which checks that a workspace has a valid instrument.
Class to represent a numeric axis of a workspace.
Definition: NumericAxis.h:29
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
Exception for errors associated with the instrument definition.
Definition: Exception.h:220
Exception for when an item is not found in a collection.
Definition: Exception.h:145
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 error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
Definition: ProgressBase.h:51
The concrete, templated class for properties.
virtual std::string value() const =0
Returns the value of the property as a string.
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
This is helper algorithm used to preprocess detector's positions namely to perform generic part of th...
bool isDetInfoLost(const Mantid::API::MatrixWorkspace_const_sptr &inWS2D) const
function checks if source workspace still has information about detectors.
void updateMasksState(const API::MatrixWorkspace_const_sptr &inputWS, DataObjects::TableWorkspace_sptr &targWS)
Method updates the column, which describes if current detector/spectra is masked It is used if one tr...
void buildFakeDetectorsPositions(const API::MatrixWorkspace_const_sptr &inputWS, DataObjects::TableWorkspace_sptr &targWS)
method calculates fake detectors positions in the situation when real detector information has been l...
std::shared_ptr< DataObjects::TableWorkspace > createTableWorkspace(const API::MatrixWorkspace_const_sptr &inputWS)
helper method to create resulting table workspace
bool m_getIsMasked
the variable specifies if one needs to return the state of detector mask e.g if the detector is maske...
double getEi(const API::MatrixWorkspace_const_sptr &inputWS) const
Method returns the efixed or Ei value stored in properties of the input workspace.
bool m_getEFixed
the variable specifies if one needs to calculate efixed for detectors (make sense for indirect instru...
void init() override
Initialize the algorithm's properties.
void processDetectorsPositions(const API::MatrixWorkspace_const_sptr &inputWS, DataObjects::TableWorkspace_sptr &targWS)
method does preliminary calculations of the detectors positions to convert results into k-dE space ; ...
void exec() override
Virtual method - must be overridden by concrete algorithm.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< TableWorkspace > TableWorkspace_sptr
shared pointer to Mantid::DataObjects::TableWorkspace
std::shared_ptr< Parameter > Parameter_sptr
Typedef for the shared pointer.
Definition: Parameter.h:195
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
Helper class which provides the Collimation Length for SANS instruments.
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
Generate a tableworkspace to store the calibration results.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54