Mantid
Loading...
Searching...
No Matches
VesuvioL1ThetaResolution.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
14#include "MantidAPI/TextAxis.h"
20#include "MantidKernel/Unit.h"
21
22#include <boost/algorithm/string/split.hpp>
23#include <boost/algorithm/string/trim.hpp>
24#include <memory>
25
26#include <fstream>
27#include <random>
28
29namespace Mantid::Algorithms {
30
31using namespace Mantid::Kernel;
32using namespace Mantid::API;
33using namespace Mantid::Geometry;
34
35namespace {
36Mantid::Kernel::Logger g_log("VesuvioL1ThetaResolution");
37}
38
39// Register the algorithm into the AlgorithmFactory
40DECLARE_ALGORITHM(VesuvioL1ThetaResolution)
41
42
43const std::string VesuvioL1ThetaResolution::name() const { return "VesuvioL1ThetaResolution"; }
44
46int VesuvioL1ThetaResolution::version() const { return 1; }
47
49const std::string VesuvioL1ThetaResolution::category() const { return "CorrectionFunctions\\SpecialCorrections"; }
50
52const std::string VesuvioL1ThetaResolution::summary() const { return "Calculates resolution of l1 and theta"; }
53
54//----------------------------------------------------------------------------------------------
58 auto positiveInt = std::make_shared<Kernel::BoundedValidator<int>>();
59 positiveInt->setLower(1);
60 auto positiveDouble = std::make_shared<Kernel::BoundedValidator<double>>();
61 positiveDouble->setLower(DBL_EPSILON);
62
63 const std::vector<std::string> exts{".par", ".dat"};
65 std::make_unique<FileProperty>("PARFile", "", FileProperty::FileAction::OptionalLoad, exts, Direction::Input),
66 "PAR file containing calibrated detector positions.");
67
68 declareProperty("SampleWidth", 3.0, positiveDouble, "With of sample in cm.");
69
70 declareProperty("SpectrumMin", 3, "Index of minimum spectrum");
71 declareProperty("SpectrumMax", 198, "Index of maximum spectrum");
72
73 declareProperty("NumEvents", 1000000, positiveInt, "Number of scattering events");
74 declareProperty("Seed", 123456789, positiveInt, "Seed for random number generator");
75
76 declareProperty("L1BinWidth", 0.05, positiveDouble, "Bin width for L1 distribution.");
77 declareProperty("ThetaBinWidth", 0.05, positiveDouble, "Bin width for theta distribution.");
78
79 declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "", Direction::Output),
80 "Output workspace containing mean and standard deviation of resolution "
81 "per detector.");
82
84 std::make_unique<WorkspaceProperty<>>("L1Distribution", "", Direction::Output, PropertyMode::Optional),
85 "Distribution of lengths of the final flight path.");
86
88 std::make_unique<WorkspaceProperty<>>("ThetaDistribution", "", Direction::Output, PropertyMode::Optional),
89 "Distribution of scattering angles.");
90}
91
92//----------------------------------------------------------------------------------------------
96 // Load the instrument workspace
98
99 const std::string l1DistributionWsName = getPropertyValue("L1Distribution");
100 const std::string thetaDistributionWsName = getPropertyValue("ThetaDistribution");
101 const size_t numHist = m_instWorkspace->getNumberHistograms();
102 const int numEvents = getProperty("NumEvents");
103
104 // Create output workspace of resolution
105 m_outputWorkspace = WorkspaceFactory::Instance().create("Workspace2D", 4, numHist, numHist);
106
107 // Set vertical axis to statistic labels
108 auto specAxis = std::make_unique<TextAxis>(4);
109 specAxis->setLabel(0, "l1_Mean");
110 specAxis->setLabel(1, "l1_StdDev");
111 specAxis->setLabel(2, "theta_Mean");
112 specAxis->setLabel(3, "theta_StdDev");
113 m_outputWorkspace->replaceAxis(1, std::move(specAxis));
114
115 // Set X axis to spectrum numbers
116 m_outputWorkspace->getAxis(0)->setUnit("Label");
117 auto xAxis = std::dynamic_pointer_cast<Units::Label>(m_outputWorkspace->getAxis(0)->unit());
118 if (xAxis)
119 xAxis->setLabel("Spectrum Number");
120
121 // Create output workspaces for distributions if required
122 if (!l1DistributionWsName.empty()) {
124
125 // Set Y axis
126 m_l1DistributionWs->setYUnitLabel("Events");
127
128 // Set X axis
129 auto distributionXAxis = m_l1DistributionWs->getAxis(0);
130 distributionXAxis->setUnit("Label");
131 auto labelUnit = std::dynamic_pointer_cast<Units::Label>(distributionXAxis->unit());
132 if (labelUnit)
133 labelUnit->setLabel("l1");
134 }
135
136 if (!thetaDistributionWsName.empty()) {
138
139 // Set Y axis
140 m_thetaDistributionWs->setYUnitLabel("Events");
141
142 // Set X axis
143 auto distributionXAxis = m_thetaDistributionWs->getAxis(0);
144 distributionXAxis->setUnit("Label");
145 auto labelUnit = std::dynamic_pointer_cast<Units::Label>(distributionXAxis->unit());
146 if (labelUnit)
147 labelUnit->setLabel("theta");
148 }
149
150 // Set up progress reporting
151 Progress prog(this, 0.0, 1.0, numHist);
152 const int seed(getProperty("Seed"));
153 std::mt19937 randEngine(static_cast<std::mt19937::result_type>(seed));
154 std::uniform_real_distribution<> flatDistrib(0.0, 1.0);
155 std::function<double()> flatVariateGen([&randEngine, &flatDistrib]() { return flatDistrib(randEngine); });
156
157 const auto &spectrumInfo = m_instWorkspace->spectrumInfo();
158 // Loop for all detectors
159 for (size_t i = 0; i < numHist; i++) {
160 std::vector<double> l1;
161 std::vector<double> theta;
162 const auto &det = spectrumInfo.detector(i);
163
164 // Report progress
165 std::stringstream report;
166 report << "Detector " << det.getID();
167 prog.report(report.str());
168 g_log.information() << "Detector ID " << det.getID() << '\n';
169
170 // Do simulation
171 calculateDetector(det, flatVariateGen, l1, theta);
172
173 // Calculate statistics for L1 and theta
174 Statistics l1Stats = getStatistics(l1);
175 Statistics thetaStats = getStatistics(theta);
176
177 g_log.information() << "l0: mean=" << l1Stats.mean << ", std.dev.=" << l1Stats.standard_deviation
178 << "\ntheta: mean=" << thetaStats.mean << ", std.dev.=" << thetaStats.standard_deviation
179 << '\n';
180
181 // Set values in output workspace
182 const int specNo = m_instWorkspace->getSpectrum(i).getSpectrumNo();
183 m_outputWorkspace->mutableX(0)[i] = specNo;
184 m_outputWorkspace->mutableX(1)[i] = specNo;
185 m_outputWorkspace->mutableX(2)[i] = specNo;
186 m_outputWorkspace->mutableX(3)[i] = specNo;
187 m_outputWorkspace->mutableY(0)[i] = l1Stats.mean;
188 m_outputWorkspace->mutableY(1)[i] = l1Stats.standard_deviation;
189 m_outputWorkspace->mutableY(2)[i] = thetaStats.mean;
190 m_outputWorkspace->mutableY(3)[i] = thetaStats.standard_deviation;
191
192 // Process data for L1 distribution
193 if (m_l1DistributionWs) {
194 auto &x = m_l1DistributionWs->mutableX(i);
195
196 std::sort(l1.begin(), l1.end());
197 std::copy(l1.begin(), l1.end(), x.begin());
198
199 m_l1DistributionWs->mutableY(i) = 1.0;
200
201 auto &spec = m_l1DistributionWs->getSpectrum(i);
202 spec.setSpectrumNo(specNo);
203 spec.addDetectorID(det.getID());
204 }
205
206 // Process data for theta distribution
208 auto &x = m_thetaDistributionWs->mutableX(i);
209
210 std::sort(theta.begin(), theta.end());
211 std::copy(theta.begin(), theta.end(), x.begin());
212
213 m_thetaDistributionWs->mutableY(i) = 1.0;
214
215 auto &spec = m_thetaDistributionWs->getSpectrum(i);
216 spec.setSpectrumNo(specNo);
217 spec.addDetectorID(det.getID());
218 }
219 }
220
221 // Process the L1 distribution workspace
222 if (m_l1DistributionWs) {
223 const double binWidth = getProperty("L1BinWidth");
224 setProperty("L1Distribution", processDistribution(m_l1DistributionWs, binWidth));
225 }
226
227 // Process the theta distribution workspace
229 const double binWidth = getProperty("ThetaBinWidth");
230 setProperty("ThetaDistribution", processDistribution(m_thetaDistributionWs, binWidth));
231 }
232
233 setProperty("OutputWorkspace", m_outputWorkspace);
234}
235
236//----------------------------------------------------------------------------------------------
240 // Get the filename for the VESUVIO IDF
241 const std::string vesuvioIPF = InstrumentFileFinder::getInstrumentFilename("VESUVIO");
242
243 // Load an empty VESUVIO instrument workspace
244 auto loadInst = AlgorithmManager::Instance().create("LoadEmptyInstrument");
245 loadInst->initialize();
246 loadInst->setChild(true);
247 loadInst->setLogging(false);
248 loadInst->setProperty("OutputWorkspace", "__evs");
249 loadInst->setProperty("Filename", vesuvioIPF);
250 loadInst->execute();
251 m_instWorkspace = loadInst->getProperty("OutputWorkspace");
252
253 // Load the PAR file if provided
254 const std::string parFilename = getPropertyValue("PARFile");
255 if (!parFilename.empty()) {
256 g_log.information() << "Loading PAR file: " << parFilename << '\n';
257
258 // Get header format
259 std::map<size_t, std::string> headerFormats;
260 headerFormats[5] = "spectrum,theta,t0,-,R";
261 headerFormats[6] = "spectrum,-,theta,t0,-,R";
262
263 std::ifstream parFile(parFilename);
264 if (!parFile) {
265 throw std::runtime_error("Cannot open PAR file");
266 }
267 std::string header;
268 getline(parFile, header);
269 g_log.debug() << "PAR file header: " << header << '\n';
270 boost::trim(header);
271 std::vector<std::string> headers;
272 boost::split(headers, header, boost::is_any_of("\t "), boost::token_compress_on);
273 size_t numCols = headers.size();
274 g_log.debug() << "PAR file columns: " << numCols << '\n';
275
276 std::string headerFormat = headerFormats[numCols];
277 if (headerFormat.empty()) {
278 std::stringstream error;
279 error << "Unrecognised PAR file header. Number of colums: " << numCols << " (expected either 5 or 6.";
280 throw std::runtime_error(error.str());
281 }
282 g_log.debug() << "PAR file header format: " << headerFormat << '\n';
283
284 // Update instrument
285 auto updateInst = AlgorithmManager::Instance().create("UpdateInstrumentFromFile");
286 updateInst->initialize();
287 updateInst->setChild(true);
288 updateInst->setLogging(false);
289 updateInst->setProperty("Workspace", m_instWorkspace);
290 updateInst->setProperty("Filename", parFilename);
291 updateInst->setProperty("MoveMonitors", false);
292 updateInst->setProperty("IgnorePhi", true);
293 updateInst->setProperty("AsciiHeader", headerFormat);
294 updateInst->execute();
295 m_instWorkspace = updateInst->getProperty("Workspace");
296 }
297
298 const int specIdxMin = static_cast<int>(m_instWorkspace->getIndexFromSpectrumNumber(getProperty("SpectrumMin")));
299 const int specIdxMax = static_cast<int>(m_instWorkspace->getIndexFromSpectrumNumber(getProperty("SpectrumMax")));
300
301 // Crop the workspace to just the detectors we are interested in
302 auto crop = AlgorithmManager::Instance().create("CropWorkspace");
303 crop->initialize();
304 crop->setChild(true);
305 crop->setLogging(false);
306 crop->setProperty("InputWorkspace", m_instWorkspace);
307 crop->setProperty("OutputWorkspace", "__evs");
308 crop->setProperty("StartWorkspaceIndex", specIdxMin);
309 crop->setProperty("EndWorkspaceIndex", specIdxMax);
310 crop->execute();
311 m_instWorkspace = crop->getProperty("OutputWorkspace");
312
313 m_sample = m_instWorkspace->getInstrument()->getSample();
314}
315
316//----------------------------------------------------------------------------------------------
320 std::function<double()> &flatRandomVariateGen,
321 std::vector<double> &l1Values, std::vector<double> &thetaValues) {
322 const int numEvents = getProperty("NumEvents");
323 l1Values.reserve(numEvents);
324 thetaValues.reserve(numEvents);
325
326 double sampleWidth = getProperty("SampleWidth");
327 // If the sample is large fix the width to the approximate beam width
328 if (sampleWidth > 4.0)
329 sampleWidth = 4.0;
330
331 // Get detector dimensions
332 Geometry::IObject_const_sptr pixelShape = detector.shape();
333 if (!pixelShape || !pixelShape->hasValidShape()) {
334 throw std::invalid_argument("Detector pixel has no defined shape!");
335 }
336 Geometry::BoundingBox detBounds = pixelShape->getBoundingBox();
337 V3D detBoxWidth = detBounds.width();
338 const double detWidth = detBoxWidth.X() * 100;
339 const double detHeight = detBoxWidth.Y() * 100;
340
341 g_log.debug() << "detWidth=" << detWidth << "\ndetHeight=" << detHeight << '\n';
342
343 // Scattering angle in rad
344 const double theta = m_instWorkspace->detectorTwoTheta(detector);
345 if (theta == 0.0)
346 return;
347
348 // Final flight path in cm
349 const double l1av = detector.getDistance(*m_sample) * 100.0;
350
351 const double x0 = l1av * sin(theta);
352 const double y0 = l1av * cos(theta);
353
354 // Get as many events as defined by NumEvents
355 // This loop is not iteration limited but highly unlikely to ever become
356 // infinate
357 while (l1Values.size() < static_cast<size_t>(numEvents)) {
358 const double xs = -sampleWidth / 2 + sampleWidth * flatRandomVariateGen();
359 const double ys = 0.0;
360 const double zs = -sampleWidth / 2 + sampleWidth * flatRandomVariateGen();
361 const double rs = sqrt(pow(xs, 2) + pow(zs, 2));
362
363 if (rs <= sampleWidth / 2) {
364 const double a = -detWidth / 2 + detWidth * flatRandomVariateGen();
365 const double xd = x0 - a * cos(theta);
366 const double yd = y0 + a * sin(theta);
367 const double zd = -detHeight / 2 + detHeight * flatRandomVariateGen();
368
369 const double l1 = sqrt(pow(xd - xs, 2) + pow(yd - ys, 2) + pow(zd - zs, 2));
370 double angle = acos(yd / l1);
371
372 if (xd < 0.0)
373 angle *= -1;
374
375 // Convert angle to degrees
376 angle *= 180.0 / M_PI;
377
378 l1Values.emplace_back(l1);
379 thetaValues.emplace_back(angle);
380 }
381
383 }
384}
385
386//----------------------------------------------------------------------------------------------
390 const size_t numHist = ws->getNumberHistograms();
391
392 double xMin(DBL_MAX);
393 double xMax(DBL_MIN);
394 for (size_t i = 0; i < numHist; i++) {
395 auto &x = ws->x(i);
396 xMin = std::min(xMin, x.front());
397 xMax = std::max(xMax, x.back());
398 }
399
400 std::stringstream binParams;
401 binParams << xMin << "," << binWidth << "," << xMax;
402
403 auto rebin = AlgorithmManager::Instance().create("Rebin");
404 rebin->initialize();
405 rebin->setChild(true);
406 rebin->setLogging(false);
407 rebin->setProperty("InputWorkspace", ws);
408 rebin->setProperty("OutputWorkspace", "__rebin");
409 rebin->setProperty("Params", binParams.str());
410 rebin->execute();
411 ws = rebin->getProperty("OutputWorkspace");
412
413 for (size_t i = 0; i < numHist; i++) {
414 auto &y = ws->y(i);
415 auto &e = ws->mutableE(i);
416 std::transform(y.begin(), y.end(), e.begin(), [](double x) { return sqrt(x); });
417 }
418
419 return ws;
420}
421
422} // namespace Mantid::Algorithms
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
double error
Definition: IndexPeaks.cpp:133
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
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
Definition: Algorithm.cpp:2026
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
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
Definition: Algorithm.cpp:1687
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
static std::string getInstrumentFilename(const std::string &instrumentName, const std::string &date="")
Get the IDF using the instrument name and date.
Helper class for reporting progress from algorithms.
Definition: Progress.h:25
A property class for workspaces.
void exec() override
Execute the algorithm.
void init() override
Initialize the algorithm's properties.
Mantid::API::MatrixWorkspace_sptr m_l1DistributionWs
int version() const override
Algorithm's version for identification.
void loadInstrument()
Loads the instrument into a workspace.
const std::string category() const override
Algorithm's category for identification.
void calculateDetector(const Mantid::Geometry::IDetector &detector, std::function< double()> &flatRandomVariateGen, std::vector< double > &l1Values, std::vector< double > &thetaValues)
Loads the instrument into a workspace.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
Mantid::API::MatrixWorkspace_sptr m_thetaDistributionWs
Mantid::API::MatrixWorkspace_sptr m_instWorkspace
Mantid::API::MatrixWorkspace_sptr processDistribution(Mantid::API::MatrixWorkspace_sptr ws, const double binWidth)
Rebins the distributions and sets error values.
Mantid::API::MatrixWorkspace_sptr m_outputWorkspace
Mantid::Geometry::IComponent_const_sptr m_sample
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
Definition: BoundingBox.h:34
Kernel::V3D width() const
Returns the width of the box.
Definition: BoundingBox.h:98
Interface class for detector objects.
Definition: IDetector.h:43
double getDistance(const IComponent &comp) const override=0
Get the distance of this detector object from another Component.
virtual const std::shared_ptr< const IObject > shape() const =0
Returns the shape of the Object.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
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
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
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
Kernel::Logger g_log("ExperimentInfo")
static logger object
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::size_t numEvents(::NeXus::File &file, bool &hasTotalCounts, bool &oldNeXusFileNames, const std::string &prefix, const NexusHDF5Descriptor &descriptor)
Get the number of events in the currently opened group.
std::shared_ptr< const IObject > IObject_const_sptr
Typdef for a shared pointer to a const object.
Definition: IObject.h:94
void MANTID_KERNEL_DLL rebin(const std::vector< double > &xold, const std::vector< double > &yold, const std::vector< double > &eold, const std::vector< double > &xnew, std::vector< double > &ynew, std::vector< double > &enew, bool distribution, bool addition=false)
Rebins data according to a new output X array.
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
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54
Simple struct to store statistics.
Definition: Statistics.h:25
double mean
Mean value.
Definition: Statistics.h:31
double standard_deviation
standard_deviation of the values
Definition: Statistics.h:35