Mantid
Loading...
Searching...
No Matches
Bin2DPowderDiffraction.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 +
24#include <algorithm>
25#include <fstream>
26#include <stdexcept>
27
28namespace Mantid::Algorithms {
29
30using namespace Kernel;
31using namespace API;
32using namespace Geometry;
33using namespace DataObjects;
34using namespace Mantid::HistogramData;
35
36// Register the algorithm into the AlgorithmFactory
37DECLARE_ALGORITHM(Bin2DPowderDiffraction)
38
39//----------------------------------------------------------------------------------------------
40
41
42const std::string Bin2DPowderDiffraction::name() const { return "Bin2DPowderDiffraction"; }
43
45int Bin2DPowderDiffraction::version() const { return 1; }
46
48const std::string Bin2DPowderDiffraction::category() const { return "Diffraction\\Focussing"; }
49
51const std::string Bin2DPowderDiffraction::summary() const {
52 return "Bins TOF powder diffraction event data in 2D space.";
53}
54
55//----------------------------------------------------------------------------------------------
59 auto wsValidator = std::make_shared<CompositeValidator>();
60 wsValidator->add<WorkspaceUnitValidator>("Wavelength");
61 wsValidator->add<SpectraAxisValidator>();
62 wsValidator->add<InstrumentValidator>();
63 wsValidator->add<HistogramValidator>();
64
66 std::make_unique<WorkspaceProperty<EventWorkspace>>("InputWorkspace", "", Direction::Input, wsValidator),
67 "An input EventWorkspace must be a Histogram workspace, not Point data. "
68 "X-axis units must be wavelength.");
69
70 declareProperty(std::make_unique<WorkspaceProperty<API::Workspace>>("OutputWorkspace", "", Direction::Output),
71 "An output workspace.");
72
73 const std::string docString = "A comma separated list of first bin boundary, width, last bin boundary. "
74 "Optionally "
75 "this can be followed by a comma and more widths and last boundary "
76 "pairs. "
77 "Negative width values indicate logarithmic binning.";
78 auto rebinValidator = std::make_shared<RebinParamsValidator>(true);
79 declareProperty(std::make_unique<ArrayProperty<double>>("dSpaceBinning", rebinValidator), docString);
80 declareProperty(std::make_unique<ArrayProperty<double>>("dPerpendicularBinning", rebinValidator), docString);
81
82 const std::vector<std::string> exts{".txt", ".dat"};
83 declareProperty(std::make_unique<FileProperty>("BinEdgesFile", "", FileProperty::OptionalLoad, exts),
84 "Optional: The ascii file containing the list of bin edges. "
85 "Either this or Axis1- and dPerpendicularBinning need to be specified.");
86
87 declareProperty(std::make_unique<PropertyWithValue<bool>>("NormalizeByBinArea", true),
88 "Normalize the binned workspace by the bin area.");
89}
90
91//----------------------------------------------------------------------------------------------
95 m_inputWS = this->getProperty("InputWorkspace");
96 m_numberOfSpectra = static_cast<int>(m_inputWS->getNumberHistograms());
97 g_log.debug() << "Number of spectra in input workspace: " << m_numberOfSpectra << "\n";
98
100
101 const bool normalizeByBinArea = this->getProperty("NormalizeByBinArea");
102 if (normalizeByBinArea)
103 normalizeToBinArea(outputWS);
104
105 setProperty("OutputWorkspace", outputWS);
106}
107//----------------------------------------------------------------------------------------------
112std::map<std::string, std::string> Bin2DPowderDiffraction::validateInputs() {
113 std::map<std::string, std::string> result;
114
115 const auto useBinFile = !getPointerToProperty("BinEdgesFile")->isDefault();
116 const auto useBinning1 = !getPointerToProperty("dSpaceBinning")->isDefault();
117 const auto useBinning2 = !getPointerToProperty("dPerpendicularBinning")->isDefault();
118 if (!useBinFile && !useBinning1 && !useBinning2) {
119 const std::string msg = "You must specify either dSpaceBinning and "
120 "dPerpendicularBinning, or a BinEdgesFile.";
121 result["dSpaceBinning"] = msg;
122 result["dPerpendicularBinning"] = msg;
123 result["BinEdgesFile"] = msg;
124 } else if (useBinFile && (useBinning1 || useBinning2)) {
125 const std::string msg = "You must specify either dSpaceBinning and "
126 "dPerpendicularBinning, or a BinEdgesFile, but not both.";
127 result["BinEdgesFile"] = msg;
128 }
129
130 return result;
131}
132//----------------------------------------------------------------------------------------------
140
142 bool binsFromFile(false);
143 size_t dPerpSize = 0;
144 size_t dSize = 0;
145 MatrixWorkspace_sptr outputWS;
146 const auto &spectrumInfo = m_inputWS->spectrumInfo();
147
148 const std::string beFileName = getProperty("BinEdgesFile");
149 if (!beFileName.empty())
150 binsFromFile = true;
151
152 const auto &oldXEdges = m_inputWS->x(0);
153 BinEdges dBins(oldXEdges.size());
154 BinEdges dPerpBins(oldXEdges.size());
155
156 auto &dPerp = dPerpBins.mutableRawData();
157 std::vector<std::vector<double>> fileXbins;
158
159 // First create the output Workspace filled with zeros
160 if (binsFromFile) {
161 dPerp.clear();
162 ReadBinsFromFile(dPerp, fileXbins);
163 dPerpSize = dPerp.size();
164 // unify xbins
165 dSize = UnifyXBins(fileXbins);
166 g_log.debug() << "Maximal size of Xbins = " << dSize;
167 outputWS = WorkspaceFactory::Instance().create(m_inputWS, dPerpSize - 1, dSize, dSize - 1);
168 g_log.debug() << "Outws has " << outputWS->getNumberHistograms() << " histograms and " << outputWS->blocksize()
169 << " bins." << std::endl;
170
171 size_t idx = 0;
172 for (const auto &Xbins : fileXbins) {
173 g_log.debug() << "Xbins size: " << Xbins.size() << std::endl;
174 BinEdges binEdges(Xbins);
175 outputWS->setBinEdges(idx, binEdges);
176 idx++;
177 }
178
179 } else {
180 static_cast<void>(createAxisFromRebinParams(getProperty("dSpaceBinning"), dBins.mutableRawData()));
181 HistogramData::BinEdges binEdges(dBins);
182 dPerpSize = createAxisFromRebinParams(getProperty("dPerpendicularBinning"), dPerp);
183 dSize = binEdges.size();
184 outputWS = WorkspaceFactory::Instance().create(m_inputWS, dPerpSize - 1, dSize, dSize - 1);
185 for (size_t idx = 0; idx < dPerpSize - 1; idx++)
186 outputWS->setBinEdges(idx, binEdges);
187 auto abscissa = std::make_unique<BinEdgeAxis>(dBins.mutableRawData());
188 outputWS->replaceAxis(0, std::move(abscissa));
189 }
190
191 outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("dSpacing");
192
193 auto verticalAxis = std::make_unique<BinEdgeAxis>(dPerp);
194 auto verticalAxisRaw = verticalAxis.get();
195 // Meta data
196 verticalAxis->unit() = UnitFactory::Instance().create("dSpacingPerpendicular");
197 verticalAxis->title() = "d_p";
198 outputWS->replaceAxis(1, std::move(verticalAxis));
199
200 Progress prog(this, 0.0, 1.0, m_numberOfSpectra);
201 auto numSpectra = static_cast<int64_t>(m_numberOfSpectra);
202 std::vector<std::vector<double>> newYValues(dPerpSize - 1, std::vector<double>(dSize - 1, 0.0));
203 std::vector<std::vector<double>> newEValues(dPerpSize - 1, std::vector<double>(dSize - 1, 0.0));
204
205 // fill the workspace with data
206 g_log.debug() << "newYSize = " << dPerpSize << std::endl;
207 g_log.debug() << "newXSize = " << dSize << std::endl;
208 std::vector<double> dp_vec(verticalAxisRaw->getValues());
209
211 for (int64_t snum = 0; snum < numSpectra; ++snum) {
213 if (!spectrumInfo.isMasked(snum)) {
214 double theta = 0.5 * spectrumInfo.twoTheta(snum);
215 double sin_theta = sin(theta);
216 if (sin_theta == 0) {
217 throw std::runtime_error("Spectrum " + std::to_string(snum) + " has sin(theta)=0. Cannot calculate d-Spacing!");
218 }
219 if (cos(theta) <= 0) {
220 throw std::runtime_error("Spectrum " + std::to_string(snum) +
221 " has cos(theta) <= 0. Cannot calculate d-SpacingPerpendicular!");
222 }
223 double log_cos_theta = log(cos(theta));
224 EventList &evList = m_inputWS->getSpectrum(snum);
225
226 // Switch to weighted if needed.
227 if (evList.getEventType() == TOF)
228 evList.switchTo(WEIGHTED);
229
230 std::vector<WeightedEvent> events = evList.getWeightedEvents();
231
232 for (const auto &ev : events) {
233 auto d = calcD(ev.tof(), sin_theta);
234 auto dp = calcDPerp(ev.tof(), log_cos_theta);
235 const auto lowy = std::lower_bound(dp_vec.begin(), dp_vec.end(), dp);
236 if ((lowy == dp_vec.end()) || (lowy == dp_vec.begin()))
237 continue;
238 int64_t dp_index = std::distance(dp_vec.begin(), lowy) - 1;
239 auto xs = binsFromFile ? fileXbins[dp_index] : dBins.rawData();
240 const auto lowx = std::lower_bound(xs.begin(), xs.end(), d);
241 if ((lowx == xs.end()) || lowx == xs.begin())
242 continue;
243 int64_t d_index = std::distance(xs.begin(), lowx) - 1;
244 // writing to the same vectors is not thread-safe
245 PARALLEL_CRITICAL(newValues) {
246 newYValues[dp_index][d_index] += ev.weight();
247 newEValues[dp_index][d_index] += ev.errorSquared();
248 }
249 }
250 }
251 prog.report("Binning event data...");
253 }
255 size_t idx = 0;
256 for (const auto &yVec : newYValues) {
257 outputWS->setCounts(idx, yVec);
258 idx++;
259 }
260 idx = 0;
261 for (auto &eVec : newEValues) {
262 std::transform(eVec.begin(), eVec.end(), eVec.begin(), static_cast<double (*)(double)>(sqrt));
263 outputWS->setCountStandardDeviations(idx, eVec);
264 idx++;
265 }
266 return outputWS;
267}
268
269//----------------------------------------------------------------------------------------------
276void Bin2DPowderDiffraction::ReadBinsFromFile(std::vector<double> &Ybins,
277 std::vector<std::vector<double>> &Xbins) const {
278 const std::string beFileName = getProperty("BinEdgesFile");
279 std::ifstream file(beFileName);
280 std::string line;
281 std::string::size_type n;
282 std::string::size_type sz;
283 std::vector<double> tmp;
284 int dpno = 0;
285 while (getline(file, line)) {
286 n = line.find("dp =");
287 if (n != std::string::npos) {
288 if (!tmp.empty()) {
289 Xbins.emplace_back(tmp);
290 tmp.clear();
291 }
292 double dp1 = std::stod(line.substr(4), &sz); // 4 is needed to crop 'dp='
293 double dp2 = std::stod(line.substr(sz + 4));
294 if (dpno < 1) {
295 Ybins.emplace_back(dp1);
296 Ybins.emplace_back(dp2);
297 } else {
298 Ybins.emplace_back(dp2);
299 }
300 dpno++;
301 } else if (line.find("#") == std::string::npos) {
302 std::stringstream ss(line);
303 double d;
304 while (ss >> d) {
305 tmp.emplace_back(d);
306 }
307 }
308 }
309 Xbins.emplace_back(tmp);
310 g_log.information() << "Number of Ybins: " << Ybins.size() << std::endl;
311 g_log.information() << "Number of Xbins sets: " << Xbins.size() << std::endl;
312}
313//----------------------------------------------------------------------------------------------
323size_t Bin2DPowderDiffraction::UnifyXBins(std::vector<std::vector<double>> &Xbins) const {
324 // get maximal vector size
325 size_t max_size = 0;
326 for (const auto &v : Xbins) {
327 max_size = std::max(v.size(), max_size);
328 }
329 // resize all vectors to maximum size, fill last vector element at the end
330 for (auto &v : Xbins) {
331 if (v.size() < max_size)
332 v.resize(max_size, v.back());
333 }
334 return max_size;
335}
336
338 NumericAxis *verticalAxis = dynamic_cast<NumericAxis *>(outWS->getAxis(1));
339 const std::vector<double> &yValues = verticalAxis->getValues();
340 auto nhist = outWS->getNumberHistograms();
341 g_log.debug() << "Number of hists: " << nhist << " Length of YAxis: " << verticalAxis->length() << std::endl;
342
343 for (size_t idx = 0; idx < nhist; ++idx) {
344 double factor = 1.0 / (yValues[idx + 1] - yValues[idx]);
345 // divide by the xBinWidth
346 outWS->convertToFrequencies(idx);
347 auto &freqs = outWS->mutableY(idx);
348 using std::placeholders::_1;
349 std::transform(freqs.begin(), freqs.end(), freqs.begin(), std::bind(std::multiplies<double>(), factor, _1));
350 auto &errors = outWS->mutableE(idx);
351 std::transform(errors.begin(), errors.end(), errors.begin(), std::bind(std::multiplies<double>(), factor, _1));
352 }
353}
354
355double calcD(double wavelength, double sintheta) { return wavelength * 0.5 / sintheta; }
356
357double calcDPerp(double wavelength, double logcostheta) { return sqrt(wavelength * wavelength - 2.0 * logcostheta); }
358
359} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
gsl_vector * tmp
#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.
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
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
A validator which checks that a workspace contains histogram data (the default) or point data as requ...
A validator which checks that a workspace has a valid instrument.
Class to represent a numeric axis of a workspace.
Definition: NumericAxis.h:29
virtual const std::vector< double > & getValues() const
Return a const reference to the values.
std::size_t length() const override
Get the length of the axis.
Definition: NumericAxis.h:38
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A validator which checks whether the input workspace has the Spectra number in the axis.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
void ReadBinsFromFile(std::vector< double > &Ybins, std::vector< std::vector< double > > &Xbins) const
Bin2DPowderDiffraction::ReadBinsFromFile.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
void normalizeToBinArea(const API::MatrixWorkspace_sptr &outWS)
void exec() override
Execute the algorithm.
int m_numberOfSpectra
The number of spectra in the workspace.
std::map< std::string, std::string > validateInputs() override
Cross-check properties with each other.
API::MatrixWorkspace_sptr createOutputWorkspace()
Setup the output workspace.
void init() override
Initialize the algorithm's properties.
int version() const override
Algorithm's version for identification.
size_t UnifyXBins(std::vector< std::vector< double > > &Xbins) const
Bin2DPowderDiffraction::UnifyXBins unifies size of the vectors in Xbins.
DataObjects::EventWorkspace_sptr m_inputWS
Pointer to the input event workspace.
const std::string category() const override
Algorithm's category for identification.
A class for holding :
Definition: EventList.h:56
Mantid::API::EventType getEventType() const override
Return the type of Event vector contained within.
Definition: EventList.cpp:643
void switchTo(Mantid::API::EventType newType) override
Switch the EventList to use the given EventType (TOF, WEIGHTED, or WEIGHTED_NOTIME)
Definition: EventList.cpp:649
std::vector< WeightedEvent > & getWeightedEvents()
Return the list of WeightedEvent contained.
Definition: EventList.cpp:795
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.
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
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 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...
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
double calcDPerp(double wavelength, double logcostheta)
double calcD(double wavelength, double sintheta)
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
STL namespace.
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