22#include <boost/algorithm/string/split.hpp>
23#include <boost/algorithm/string/trim.hpp>
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);
63 const std::vector<std::string> exts{
".par",
".dat"};
66 "PAR file containing calibrated detector positions.");
68 declareProperty(
"SampleWidth", 3.0, positiveDouble,
"With of sample in cm.");
73 declareProperty(
"NumEvents", 1000000, positiveInt,
"Number of scattering events");
74 declareProperty(
"Seed", 123456789, positiveInt,
"Seed for random number generator");
76 declareProperty(
"L1BinWidth", 0.05, positiveDouble,
"Bin width for L1 distribution.");
77 declareProperty(
"ThetaBinWidth", 0.05, positiveDouble,
"Bin width for theta distribution.");
80 "Output workspace containing mean and standard deviation of resolution "
85 "Distribution of lengths of the final flight path.");
89 "Distribution of scattering angles.");
100 const std::string thetaDistributionWsName =
getPropertyValue(
"ThetaDistribution");
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");
117 auto xAxis = std::dynamic_pointer_cast<Units::Label>(
m_outputWorkspace->getAxis(0)->unit());
119 xAxis->setLabel(
"Spectrum Number");
122 if (!l1DistributionWsName.empty()) {
130 distributionXAxis->setUnit(
"Label");
131 auto labelUnit = std::dynamic_pointer_cast<Units::Label>(distributionXAxis->unit());
133 labelUnit->setLabel(
"l1");
136 if (!thetaDistributionWsName.empty()) {
144 distributionXAxis->setUnit(
"Label");
145 auto labelUnit = std::dynamic_pointer_cast<Units::Label>(distributionXAxis->unit());
147 labelUnit->setLabel(
"theta");
151 Progress prog(
this, 0.0, 1.0, numHist);
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); });
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);
165 std::stringstream report;
166 report <<
"Detector " << det.getID();
167 prog.
report(report.str());
196 std::sort(l1.begin(), l1.end());
197 std::copy(l1.begin(), l1.end(),
x.begin());
202 spec.setSpectrumNo(specNo);
203 spec.addDetectorID(det.getID());
210 std::sort(theta.begin(), theta.end());
211 std::copy(theta.begin(), theta.end(),
x.begin());
216 spec.setSpectrumNo(specNo);
217 spec.addDetectorID(det.getID());
229 const double binWidth =
getProperty(
"ThetaBinWidth");
245 loadInst->initialize();
246 loadInst->setChild(
true);
247 loadInst->setLogging(
false);
248 loadInst->setProperty(
"OutputWorkspace",
"__evs");
249 loadInst->setProperty(
"Filename", vesuvioIPF);
255 if (!parFilename.empty()) {
259 std::map<size_t, std::string> headerFormats;
260 headerFormats[5] =
"spectrum,theta,t0,-,R";
261 headerFormats[6] =
"spectrum,-,theta,t0,-,R";
263 std::ifstream parFile(parFilename);
265 throw std::runtime_error(
"Cannot open PAR file");
268 getline(parFile, header);
269 g_log.
debug() <<
"PAR file header: " << header <<
'\n';
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';
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());
282 g_log.
debug() <<
"PAR file header format: " << headerFormat <<
'\n';
286 updateInst->initialize();
287 updateInst->setChild(
true);
288 updateInst->setLogging(
false);
290 updateInst->setProperty(
"Filename", parFilename);
291 updateInst->setProperty(
"MoveMonitors",
false);
292 updateInst->setProperty(
"IgnorePhi",
true);
293 updateInst->setProperty(
"AsciiHeader", headerFormat);
294 updateInst->execute();
304 crop->setChild(
true);
305 crop->setLogging(
false);
307 crop->setProperty(
"OutputWorkspace",
"__evs");
308 crop->setProperty(
"StartWorkspaceIndex", specIdxMin);
309 crop->setProperty(
"EndWorkspaceIndex", specIdxMax);
320 std::function<
double()> &flatRandomVariateGen,
321 std::vector<double> &l1Values, std::vector<double> &thetaValues) {
328 if (sampleWidth > 4.0)
333 if (!pixelShape || !pixelShape->hasValidShape()) {
334 throw std::invalid_argument(
"Detector pixel has no defined shape!");
337 V3D detBoxWidth = detBounds.
width();
338 const double detWidth = detBoxWidth.
X() * 100;
339 const double detHeight = detBoxWidth.
Y() * 100;
341 g_log.
debug() <<
"detWidth=" << detWidth <<
"\ndetHeight=" << detHeight <<
'\n';
351 const double x0 = l1av * sin(theta);
352 const double y0 = l1av * cos(theta);
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));
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();
369 const double l1 = sqrt(pow(xd - xs, 2) + pow(yd - ys, 2) + pow(zd - zs, 2));
370 double angle = acos(yd / l1);
376 angle *= 180.0 / M_PI;
378 l1Values.emplace_back(l1);
379 thetaValues.emplace_back(angle);
390 const size_t numHist = ws->getNumberHistograms();
392 double xMin(DBL_MAX);
393 double xMax(DBL_MIN);
394 for (
size_t i = 0; i < numHist; i++) {
396 xMin = std::min(xMin,
x.front());
397 xMax = std::max(xMax,
x.back());
400 std::stringstream binParams;
401 binParams << xMin <<
"," << binWidth <<
"," << xMax;
405 rebin->setChild(
true);
406 rebin->setLogging(
false);
407 rebin->setProperty(
"InputWorkspace", ws);
408 rebin->setProperty(
"OutputWorkspace",
"__rebin");
409 rebin->setProperty(
"Params", binParams.str());
411 ws =
rebin->getProperty(
"OutputWorkspace");
413 for (
size_t i = 0; i < numHist; i++) {
415 auto &e = ws->mutableE(i);
416 std::transform(
y.begin(),
y.end(), e.begin(), [](
double x) { return sqrt(x); });
#define DECLARE_ALGORITHM(classname)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
void interruption_point()
This is called during long-running operations, and check if the algorithm has requested that it be ca...
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
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.
A property class for workspaces.
VesuvioL1ThetaResolution.
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.
Kernel::V3D width() const
Returns the width of the box.
Interface class for detector objects.
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.
void debug(const std::string &msg)
Logs at debug level.
void information(const std::string &msg)
Logs at information level.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
constexpr double X() const noexcept
Get x.
constexpr double Y() const noexcept
Get y.
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.
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.
@ Input
An input workspace.
@ Output
An output workspace.
Simple struct to store statistics.
double standard_deviation
standard_deviation of the values