Mantid
Loading...
Searching...
No Matches
CreateSimulationWorkspace.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/Axis.h"
11#include "MantidAPI/Run.h"
14
23
24#include "LoadRaw/isisraw2.h"
25// clang-format off
26#include <nexus/NeXusFile.hpp>
27#include <nexus/NeXusException.hpp>
28// clang-format on
29
30#include <Poco/File.h>
31
32#include <boost/algorithm/string/predicate.hpp>
33
34namespace {
35
36struct StartAndEndTime {
37 Mantid::Types::Core::DateAndTime startTime;
38 Mantid::Types::Core::DateAndTime endTime;
39};
40
41StartAndEndTime getStartAndEndTimesFromRawFile(const std::string &filename) {
42 FILE *rawFile = fopen(filename.c_str(), "rb");
43 if (!rawFile)
44 throw std::runtime_error("Cannot open RAW file for reading: " + filename);
45
46 ISISRAW2 isisRaw;
47 const bool fromFile(true), readData(false);
48 isisRaw.ioRAW(rawFile, fromFile, readData);
49
50 StartAndEndTime startAndEndTime;
51 startAndEndTime.startTime = Mantid::DataHandling::LoadRawHelper::extractStartTime(isisRaw);
52 startAndEndTime.endTime = Mantid::DataHandling::LoadRawHelper::extractEndTime(isisRaw);
53
54 fclose(rawFile);
55 return startAndEndTime;
56}
57
58StartAndEndTime getStartAndEndTimesFromNexusFile(const std::string &filename,
59 const Mantid::Types::Core::DateAndTime &startTimeDefault,
60 const Mantid::Types::Core::DateAndTime &endTimeDefault) {
61 StartAndEndTime startAndEndTime;
62 try {
63 startAndEndTime.startTime = Mantid::DataHandling::extractStartTime(filename);
64 startAndEndTime.endTime = Mantid::DataHandling::extractEndTime(filename);
65 } catch (...) {
66 startAndEndTime.startTime = startTimeDefault;
67 startAndEndTime.endTime = endTimeDefault;
68 }
69
70 return startAndEndTime;
71}
72} // namespace
73
74namespace Mantid::DataHandling {
75
76// Register the algorithm into the AlgorithmFactory
77DECLARE_ALGORITHM(CreateSimulationWorkspace)
78
79using namespace API;
80using namespace HistogramData;
81
82//----------------------------------------------------------------------------------------------
84const std::string CreateSimulationWorkspace::name() const { return "CreateSimulationWorkspace"; }
85
87int CreateSimulationWorkspace::version() const { return 1; }
88
90const std::string CreateSimulationWorkspace::category() const { return "Utility\\Workspaces"; }
91
92//----------------------------------------------------------------------------------------------
93
94//----------------------------------------------------------------------------------------------
98 using namespace Kernel;
99
100 declareProperty("Instrument", "", std::make_shared<MandatoryValidator<std::string>>(),
101 "An instrument name or filename ( a full path or string "
102 "containing an xml extension).",
104
106 std::make_unique<ArrayProperty<double>>("BinParams", std::make_shared<RebinParamsValidator>(), Direction::Input),
107 "A comma separated list of first bin boundary, width, last "
108 "bin boundary. See Rebin for more details");
109
110 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output), "The new workspace");
111
112 auto knownUnits = UnitFactory::Instance().getKeys();
113 declareProperty("UnitX", "DeltaE", std::make_shared<ListValidator<std::string>>(knownUnits),
114 "The unit to assign to the X axis", Direction::Input);
115
117 std::make_unique<FileProperty>("DetectorTableFilename", "", FileProperty::OptionalLoad, "", Direction::Input),
118 "An optional filename (currently RAW or ISIS NeXus) that "
119 "contains UDET & SPEC tables to access hardware grouping");
120
121 declareProperty("SetErrors", false, "Whether to set histogram bin errors to sqrt of intensity.");
122}
123
124//----------------------------------------------------------------------------------------------
130
131 setProperty("OutputWorkspace", m_outputWS);
132}
133
134//-----------------------------------------------------------------------------------------------
135// Private methods
136//-----------------------------------------------------------------------------------------------
142 const bool enableLogging(false);
143 auto loadInstrument = createChildAlgorithm("LoadInstrument", 0.0, 0.5, enableLogging);
144 MatrixWorkspace_sptr tempWS = WorkspaceFactory::Instance().create("Workspace2D", 1, 1, 1);
145
146 // We need to set the correct start date for this workspace
147 // else we might be pulling an inadequate IDF
148 setStartDate(tempWS);
149
150 loadInstrument->setProperty("Workspace", tempWS);
151 const std::string instrProp = getProperty("Instrument");
152 if (boost::algorithm::ends_with(instrProp, ".xml")) {
153 loadInstrument->setPropertyValue("Filename", instrProp);
154 } else {
155 loadInstrument->setPropertyValue("InstrumentName", instrProp);
156 }
157 loadInstrument->setProperty("RewriteSpectraMap", Kernel::OptionalBool(true));
158 loadInstrument->executeAsChildAlg();
159 tempWS = loadInstrument->getProperty("Workspace");
160
161 m_instrument = tempWS->getInstrument();
162}
163
168 const size_t nhistograms = createDetectorMapping();
169 const auto binBoundaries = createBinBoundaries();
170 const size_t xlength = binBoundaries.size();
171 const size_t ylength = xlength - 1;
172 const bool setError = getProperty("SetErrors");
173
174 m_outputWS = WorkspaceFactory::Instance().create("Workspace2D", nhistograms, xlength, ylength);
175 m_outputWS->setInstrument(m_instrument);
176 m_outputWS->populateInstrumentParameters();
177
178 m_outputWS->getAxis(0)->setUnit(getProperty("UnitX"));
179 m_outputWS->setYUnit("SpectraNumber");
180
181 m_progress = std::make_shared<Progress>(this, 0.5, 0.75, nhistograms);
182
184 for (int64_t i = 0; i < static_cast<int64_t>(nhistograms); ++i) {
185 m_outputWS->setBinEdges(i, binBoundaries);
186 m_outputWS->mutableY(i) = 1.0;
187 if (setError) {
188 m_outputWS->mutableE(i).assign(m_outputWS->getNumberBins(i), sqrt(1.0));
189 }
190
191 m_progress->report("Setting X values");
192 }
194
195 // Update the instrument from the file if necessary
196 const std::string detTableFile = getProperty("DetectorTableFilename");
197 if (boost::algorithm::ends_with(detTableFile, ".raw") || boost::algorithm::ends_with(detTableFile, ".RAW") ||
198 boost::algorithm::ends_with(detTableFile, ".nxs") || boost::algorithm::ends_with(detTableFile, ".NXS")) {
199 adjustInstrument(detTableFile);
200 }
201}
202
209 const std::string detTableFile = getProperty("DetectorTableFilename");
210 if (detTableFile.empty()) {
212 } else {
213 loadMappingFromFile(detTableFile);
214 }
215 return m_detGroups.size();
216}
217
222 const std::vector<detid_t> detids = m_instrument->getDetectorIDs(true);
223 const size_t nhist = detids.size();
224
225 m_detGroups.clear();
226 for (size_t i = 0; i < nhist; ++i) {
227 std::set<detid_t> group;
228 group.insert(detids[i]);
229 m_detGroups.emplace(static_cast<specnum_t>(i + 1), group);
230 }
231}
232
237void CreateSimulationWorkspace::loadMappingFromFile(const std::string &filename) {
238 if (boost::algorithm::ends_with(filename, ".raw") || boost::algorithm::ends_with(filename, ".RAW")) {
239 loadMappingFromRAW(filename);
240 } else if (boost::algorithm::ends_with(filename, ".nxs") || boost::algorithm::ends_with(filename, ".NXS")) {
241 loadMappingFromISISNXS(filename);
242 }
243}
244
249void CreateSimulationWorkspace::loadMappingFromRAW(const std::string &filename) {
250 FILE *rawFile = fopen(filename.c_str(), "rb");
251 if (!rawFile)
252 throw std::runtime_error("Cannot open RAW file for reading: " + filename);
253
254 ISISRAW2 isisRaw;
255 const bool fromFile(true), readData(false);
256 isisRaw.ioRAW(rawFile, fromFile, readData);
257
258 int ndet = isisRaw.i_det;
259 int *specTable = isisRaw.spec;
260 int *udetTable = isisRaw.udet;
261 createGroupingsFromTables(specTable, udetTable, ndet);
262
263 fclose(rawFile);
264}
265
272void CreateSimulationWorkspace::loadMappingFromISISNXS(const std::string &filename) {
273 ::NeXus::File nxsFile(filename);
274 try {
275 nxsFile.openPath("/raw_data_1/isis_vms_compat");
276 } catch (::NeXus::Exception &) {
277 throw std::runtime_error("Cannot find path to isis_vms_compat. Is the file an ISIS NeXus file?");
278 }
279 using NXIntArray = std::unique_ptr<std::vector<int32_t>>;
280
281 nxsFile.openData("NDET");
282 NXIntArray ndets(nxsFile.getData<int32_t>());
283 nxsFile.closeData();
284
285 nxsFile.openData("SPEC");
286 NXIntArray specTable(nxsFile.getData<int32_t>());
287 nxsFile.closeData();
288
289 nxsFile.openData("UDET");
290 NXIntArray udetTable(nxsFile.getData<int32_t>());
291 nxsFile.closeData();
292
293 createGroupingsFromTables(specTable->data(), udetTable->data(), (*ndets)[0]);
294}
295
302void CreateSimulationWorkspace::createGroupingsFromTables(const int *specTable, const int *udetTable, int ndets) {
303 m_detGroups.clear();
304 for (int i = 0; i < ndets; ++i) {
305 int specNo = specTable[i];
306 int detID = udetTable[i];
307 if (m_instrument->isMonitor(detID))
308 continue; // Skip monitors
309
310 auto iter = m_detGroups.find(specNo);
311 if (iter != m_detGroups.end()) {
312 iter->second.insert(detID);
313 } else {
314 std::set<detid_t> group;
315 group.insert(static_cast<detid_t>(detID));
316 m_detGroups.emplace(specNo, group);
317 }
318 }
319}
320
325 const std::vector<double> rbparams = getProperty("BinParams");
326 MantidVec newBins;
327 const int numBoundaries = Mantid::Kernel::VectorHelper::createAxisFromRebinParams(rbparams, newBins);
328 if (numBoundaries <= 2) {
329 throw std::invalid_argument("Error in BinParams - Gave invalid number of bin boundaries: " +
330 std::to_string(numBoundaries));
331 }
332 return BinEdges(std::move(newBins));
333}
334
339 size_t wsIndex(0);
340 for (auto &detGroup : m_detGroups) {
341 auto &spectrum = m_outputWS->getSpectrum(wsIndex);
342 spectrum.setSpectrumNo(static_cast<specnum_t>(wsIndex + 1)); // Ensure a contiguous mapping
343 spectrum.clearDetectorIDs();
344 spectrum.addDetectorIDs(detGroup.second);
345 ++wsIndex;
346 }
347}
348
353void CreateSimulationWorkspace::adjustInstrument(const std::string &filename) {
354 // If requested update the instrument to positions in the raw file
355 const auto &pmap = m_outputWS->constInstrumentParameters();
356 Geometry::Instrument_const_sptr instrument = m_outputWS->getInstrument();
357 std::shared_ptr<Geometry::Parameter> updateDets = pmap.get(instrument->getComponentID(), "det-pos-source");
358 if (!updateDets)
359 return; // No tag, use IDF
360
361 std::string value = updateDets->value<std::string>();
362 if (value.substr(0, 8) == "datafile") {
363 auto updateInst = createChildAlgorithm("UpdateInstrumentFromFile", 0.75, 1.0);
364 updateInst->setProperty<MatrixWorkspace_sptr>("Workspace", m_outputWS);
365 updateInst->setPropertyValue("Filename", filename);
366 if (value == "datafile-ignore-phi") {
367 updateInst->setProperty("IgnorePhi", true);
368 g_log.information("Detector positions in IDF updated with positions in "
369 "the data file except for the phi values");
370 } else {
371 g_log.information("Detector positions in IDF updated with positions in the data file");
372 }
373 // We want this to throw if it fails to warn the user that the information
374 // is not correct.
375 updateInst->execute();
376 }
377}
378
385 const std::string detTableFile = getProperty("DetectorTableFilename");
386 auto hasDetTableFile = !detTableFile.empty();
387 auto &run = workspace->mutableRun();
388
389 Types::Core::DateAndTime startTime;
390 Types::Core::DateAndTime endTime;
391 try {
392 // The start and end times might not be valid, and hence can throw
393 startTime = run.startTime();
394 endTime = run.endTime();
395 } catch (std::runtime_error &) {
396 startTime = Types::Core::DateAndTime::getCurrentTime();
397 endTime = Types::Core::DateAndTime::getCurrentTime();
398 }
399
400 if (hasDetTableFile) {
401 if (boost::algorithm::ends_with(detTableFile, ".raw") || boost::algorithm::ends_with(detTableFile, ".RAW")) {
402 auto startAndEndTime = getStartAndEndTimesFromRawFile(detTableFile);
403 startTime = startAndEndTime.startTime;
404 endTime = startAndEndTime.endTime;
405 } else if (boost::algorithm::ends_with(detTableFile, ".nxs") || boost::algorithm::ends_with(detTableFile, ".NXS")) {
406 auto startAndEndTime = getStartAndEndTimesFromNexusFile(detTableFile, startTime, endTime);
407 startTime = startAndEndTime.startTime;
408 endTime = startAndEndTime.endTime;
409 }
410 }
411
412 run.setStartAndEndTime(startTime, endTime);
413}
414
415} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double value
The value of the point.
Definition: FitMW.cpp:51
IPeaksWorkspace_sptr workspace
Definition: IndexPeaks.cpp:114
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
isis raw file.
Definition: isisraw2.h:13
int ioRAW(FILE *file, bool from_file, bool read_data=true) override
Loads the headers of the file, leaves the file pointer at a specific position.
Definition: isisraw2.cpp:43
int i_det
number of detectors NDET
Definition: isisraw.h:296
int * udet
user detector number for each detector (size NDET)
Definition: isisraw.h:322
int * spec
spectrum number table (size NDET)
Definition: isisraw.h:302
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
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
A property class for workspaces.
API::MatrixWorkspace_sptr m_outputWS
Pointer to the new workspace.
void createGroupingsFromTables(const int *specTable, const int *udetTable, int ndets)
Create the grouping map from the tables.
void init() override
Initialize the algorithm's properties.
std::shared_ptr< API::Progress > m_progress
Pointer to a progress object.
std::map< specnum_t, std::set< detid_t > > m_detGroups
List of detector groupings.
void createOutputWorkspace()
Creates the output workspace.
void setStartDate(const API::MatrixWorkspace_sptr &workspace)
Set start date for dummy workspace.
const std::string name() const override
Algorithm's name for identification.
int version() const override
Algorithm's version for identification.
void loadMappingFromFile(const std::string &filename)
Load the detector mapping from a file.
void applyDetectorMapping()
Apply the created mapping to the workspace.
void createOneToOneMapping()
Create a one to one mapping from the spectrum numbers to detector IDs.
void loadMappingFromISISNXS(const std::string &filename)
Load the detector mapping from a NXS file.
Geometry::Instrument_const_sptr m_instrument
Pointer to the new instrument.
size_t createDetectorMapping()
Creates the detector grouping list.
HistogramData::BinEdges createBinBoundaries() const
Returns new Xbins.
const std::string category() const override
Algorithm's category for identification.
void adjustInstrument(const std::string &filename)
Apply any instrument adjustments from the file.
void loadMappingFromRAW(const std::string &filename)
Load the detector mapping from a RAW file.
static Types::Core::DateAndTime extractEndTime(ISISRAW &isisRaw)
Extract the end time from a raw file.
static Types::Core::DateAndTime extractStartTime(ISISRAW &isisRaw)
Extract the start time from a raw file.
Support for a property that holds an array of values.
Definition: ArrayProperty.h:28
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
ListValidator is a validator that requires the value of a property to be one of a defined list of pos...
Definition: ListValidator.h:29
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Validator to check that a property is not left empty.
OptionalBool : Tri-state bool.
Definition: OptionalBool.h:25
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
Mantid::Types::Core::DateAndTime DLLExport extractEndTime(const std::string &filename)
Gets the start time from the nexus file.
Mantid::Types::Core::DateAndTime DLLExport extractStartTime(const std::string &filename)
Extracts the start and the end time from a Nexus file.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
int MANTID_KERNEL_DLL createAxisFromRebinParams(const std::vector< double > &params, std::vector< double > &xnew, const bool resize_xnew=true, const bool full_bins_only=false, const double xMinHint=std::nan(""), const double xMaxHint=std::nan(""), const bool useReverseLogarithmic=false, const double power=-1)
Creates a new output X array given a 'standard' set of rebinning parameters.
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
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
Definition: cow_ptr.h:172
int32_t specnum_t
Typedef for a spectrum Number.
Definition: IDTypes.h:16
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54