Mantid
Loading...
Searching...
No Matches
GenerateGroupingPowder.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 +
16#include "MantidKernel/System.h"
17
18#include <Poco/DOM/AutoPtr.h>
19#include <Poco/DOM/DOMWriter.h>
20#include <Poco/DOM/Document.h>
21#include <Poco/DOM/Element.h>
22#include <Poco/DOM/Text.h>
23#include <Poco/XML/XMLWriter.h>
24
25#include <fstream>
26
27using namespace Mantid::API;
28using namespace Mantid::Kernel;
29using namespace Mantid::Geometry;
30using namespace Poco::XML;
31
32namespace Mantid::DataHandling {
33
34// Register the algorithm into the AlgorithmFactory
36
37
38const std::string GenerateGroupingPowder::name() const { return "GenerateGroupingPowder"; }
39
41int GenerateGroupingPowder::version() const { return 1; }
42
44const std::string GenerateGroupingPowder::category() const {
45 return R"(DataHandling\Grouping;Transforms\Grouping;Diffraction\Utility)";
46}
47
51 declareProperty(std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input),
52 "A workspace from which to generate the grouping.");
53 auto positiveDouble = std::make_shared<BoundedValidator<double>>();
54 positiveDouble->setLower(0.0);
55 declareProperty("AngleStep", -1.0, positiveDouble, "The angle step for grouping, in degrees.");
56 declareProperty(std::make_unique<FileProperty>("GroupingFilename", "", FileProperty::Save, ".xml"),
57 "A grouping file that will be created.");
58 declareProperty("GenerateParFile", true,
59 "If true, a par file with a corresponding name to the "
60 "grouping file will be generated.");
61}
62
66 MatrixWorkspace_const_sptr input_ws = getProperty("InputWorkspace");
67 const auto &spectrumInfo = input_ws->spectrumInfo();
68 const auto &detectorInfo = input_ws->detectorInfo();
69 const auto &detectorIDs = detectorInfo.detectorIDs();
70 if (detectorIDs.empty())
71 throw std::invalid_argument("Workspace contains no detectors.");
72
73 const double step = getProperty("AngleStep");
74 const auto numSteps = static_cast<size_t>(180. / step + 1);
75
76 std::vector<std::vector<detid_t>> groups(numSteps);
77 std::vector<double> twoThetaAverage(numSteps, 0.);
78 std::vector<double> rAverage(numSteps, 0.);
79
80 for (size_t i = 0; i < spectrumInfo.size(); ++i) {
81 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMasked(i) || spectrumInfo.isMonitor(i)) {
82 continue;
83 }
84 const auto &det = spectrumInfo.detector(i);
85 const double tt = spectrumInfo.twoTheta(i) * Geometry::rad2deg;
86 const double r = spectrumInfo.l2(i);
87 const auto where = static_cast<size_t>(tt / step);
88 twoThetaAverage[where] += tt;
89 rAverage[where] += r;
90 if (spectrumInfo.hasUniqueDetector(i)) {
91 groups[where].emplace_back(det.getID());
92 } else {
93 const auto &group = dynamic_cast<const DetectorGroup &>(det);
94 const auto ids = group.getDetectorIDs();
95 groups[where].insert(groups[where].end(), ids.begin(), ids.end());
96 }
97 }
98
99 const std::string XMLfilename = getProperty("GroupingFilename");
100
101 // XML
102 AutoPtr<Document> pDoc = new Document;
103 AutoPtr<Element> pRoot = pDoc->createElement("detector-grouping");
104 pDoc->appendChild(pRoot);
105 pRoot->setAttribute("instrument", input_ws->getInstrument()->getName());
106
107 size_t goodGroups(0);
108 for (size_t i = 0; i < numSteps; ++i) {
109 size_t gSize = groups.at(i).size();
110 if (gSize > 0) {
111 ++goodGroups;
112 std::stringstream spID, textvalue;
113 spID << i;
114 AutoPtr<Element> pChildGroup = pDoc->createElement("group");
115 pChildGroup->setAttribute("ID", spID.str());
116 pRoot->appendChild(pChildGroup);
117
118 std::copy(groups.at(i).begin(), groups.at(i).end(), std::ostream_iterator<size_t>(textvalue, ","));
119 std::string text = textvalue.str();
120 const size_t found = text.rfind(',');
121 if (found != std::string::npos) {
122 text.erase(found, 1); // erase the last comma
123 }
124
125 AutoPtr<Element> pDetid = pDoc->createElement("detids");
126 AutoPtr<Text> pText1 = pDoc->createTextNode(text);
127 pDetid->appendChild(pText1);
128 pChildGroup->appendChild(pDetid);
129 }
130 }
131 if (goodGroups == 0) {
132 throw Exception::InstrumentDefinitionError("No detectors found in scattering angles between 0 and 180 degrees");
133 }
134
135 DOMWriter writer;
136 writer.setNewLine("\n");
137 writer.setOptions(XMLWriter::PRETTY_PRINT);
138
139 std::ofstream ofs(XMLfilename.c_str());
140 if (!ofs) {
141 g_log.error("Unable to create file: " + XMLfilename);
142 throw Exception::FileError("Unable to create file: ", XMLfilename);
143 }
144 ofs << "<?xml version=\"1.0\"?>\n";
145
146 writer.writeNode(ofs, pDoc);
147 ofs.close();
148
149 // PAR file
150 const bool generatePar = getProperty("GenerateParFile");
151 if (generatePar) {
152 std::string PARfilename = XMLfilename;
153 PARfilename.replace(PARfilename.end() - 3, PARfilename.end(), "par");
154 std::ofstream outPAR_file(PARfilename.c_str());
155 if (!outPAR_file) {
156 g_log.error("Unable to create file: " + PARfilename);
157 throw Exception::FileError("Unable to create file: ", PARfilename);
158 }
159 // Write the number of detectors to the file.
160 outPAR_file << " " << goodGroups << '\n';
161
162 for (size_t i = 0; i < numSteps; ++i) {
163 const size_t gSize = groups.at(i).size();
164 if (gSize > 0) {
165 outPAR_file << std::fixed << std::setprecision(3);
166 outPAR_file.width(10);
167 outPAR_file << rAverage.at(i) / static_cast<double>(gSize);
168 outPAR_file.width(10);
169 outPAR_file << twoThetaAverage.at(i) / static_cast<double>(gSize);
170 outPAR_file.width(10);
171 outPAR_file << 0.;
172 outPAR_file.width(10);
173 outPAR_file << step * Geometry::deg2rad * rAverage.at(i) / static_cast<double>(gSize);
174 outPAR_file.width(10);
175 outPAR_file << 0.01;
176 outPAR_file.width(10);
177 outPAR_file << (groups.at(i)).at(0) << '\n';
178 }
179 }
180
181 // Close the file
182 outPAR_file.close();
183 }
184}
185
186} // namespace Mantid::DataHandling
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
@ Save
to specify a file to write to, the file may or may not exist
Definition: FileProperty.h:49
A property class for workspaces.
GenerateGroupingPowder : Generate grouping file and par file, for powder scattering.
int version() const override
Algorithm's version for identification.
void exec() override
Execute the algorithm.
const std::string category() const override
Algorithm's category for identification.
void init() override
Initialize the algorithm's properties.
Holds a collection of detectors.
Definition: DetectorGroup.h:28
std::vector< detid_t > getDetectorIDs() const
What detectors are contained in the group?
Records the filename and the description of failure.
Definition: Exception.h:98
Exception for errors associated with the instrument definition.
Definition: Exception.h:220
void error(const std::string &msg)
Logs at error level.
Definition: Logger.cpp:77
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
constexpr double deg2rad
Defines units/enum for Crystal work.
Definition: AngleUnits.h:20
constexpr double rad2deg
Radians to degrees conversion factor.
Definition: AngleUnits.h:23
STL namespace.
@ Input
An input workspace.
Definition: Property.h:53