Mantid
Loading...
Searching...
No Matches
LoadIsawSpectrum.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"
17#include "MantidKernel/Unit.h"
19#include "MantidKernel/Utils.h"
20
21#include <fstream>
22
23using namespace Mantid::Geometry;
24using namespace Mantid::DataObjects;
25using namespace Mantid::Kernel;
26using namespace Mantid::API;
27using namespace Mantid::PhysicalConstants;
28
29namespace Mantid::Crystal {
30
31// Register the algorithm into the AlgorithmFactory
32DECLARE_ALGORITHM(LoadIsawSpectrum)
33
34
36void LoadIsawSpectrum::init() {
37 declareProperty(std::make_unique<FileProperty>("SpectraFile", "", API::FileProperty::Load, ".dat"),
38 "Incident spectrum and detector efficiency correction file.");
39 declareProperty(std::make_unique<WorkspaceProperty<MatrixWorkspace>>("OutputWorkspace", "", Direction::Output),
40 "An output Workspace containing spectra for each detector bank.");
41 // 3 properties for getting the right instrument
42 getInstrument3WaysInit(this);
43}
44
49
50 // If sample not at origin, shift cached positions.
51 const V3D samplePos = inst->getSample()->getPos();
52 const V3D pos = inst->getSource()->getPos() - samplePos;
53 double l1 = pos.norm();
54
55 std::string STRING;
56 std::ifstream infile;
57 std::string spectraFile = getPropertyValue("SpectraFile");
58 infile.open(spectraFile.c_str());
59
60 size_t a = -1;
61 std::vector<std::vector<double>> spectra;
62 std::vector<std::vector<double>> time;
63 int iSpec = 0;
64 for (int wi = 0; wi < 8; wi++)
65 getline(infile, STRING); // Saves the line in STRING.
66 while (!infile.eof()) // To get you all the lines.
67 {
68 time.resize(a + 1);
69 spectra.resize(a + 1);
70 getline(infile, STRING); // Saves the line in STRING.
71 if (infile.eof())
72 break;
73 std::stringstream ss(STRING);
74 if (STRING.find("Bank") == std::string::npos) {
75 double time0, spectra0;
76 ss >> time0 >> spectra0;
77 time[a].emplace_back(time0);
78 spectra[a].emplace_back(spectra0);
79
80 } else {
81 a++;
82 }
83 }
84 infile.close();
85 // Build a list of Rectangular Detectors
86 std::vector<std::shared_ptr<RectangularDetector>> detList;
87 for (int i = 0; i < inst->nelements(); i++) {
88 std::shared_ptr<RectangularDetector> det;
89 std::shared_ptr<ICompAssembly> assem;
90 std::shared_ptr<ICompAssembly> assem2;
91
92 det = std::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
93 if (det) {
94 detList.emplace_back(det);
95 } else {
96 // Also, look in the first sub-level for RectangularDetectors (e.g. PG3).
97 // We are not doing a full recursive search since that will be very long
98 // for lots of pixels.
99 assem = std::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
100 if (assem) {
101 for (int j = 0; j < assem->nelements(); j++) {
102 det = std::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
103 if (det) {
104 detList.emplace_back(det);
105 } else {
106 // Also, look in the second sub-level for RectangularDetectors (e.g.
107 // PG3).
108 // We are not doing a full recursive search since that will be very
109 // long for lots of pixels.
110 assem2 = std::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
111 if (assem2) {
112 for (int k = 0; k < assem2->nelements(); k++) {
113 det = std::dynamic_pointer_cast<RectangularDetector>((*assem2)[k]);
114 if (det) {
115 detList.emplace_back(det);
116 }
117 }
118 }
119 }
120 }
121 }
122 }
123 }
124
125 if (spectra.size() < 1)
126 throw std::runtime_error("The number of spectra in the loaded file is zero.");
127
128 MatrixWorkspace_sptr outWS = std::dynamic_pointer_cast<MatrixWorkspace>(
129 API::WorkspaceFactory::Instance().create("Workspace2D", spectra.size(), spectra[0].size(), spectra[0].size()));
130 outWS->setInstrument(inst);
131 outWS->getAxis(0)->setUnit("TOF");
132 outWS->setYUnit("Counts");
133 outWS->setDistribution(true);
134 outWS->rebuildSpectraMapping(false);
135
136 // Go through each point at this run / bank
137 for (size_t i = 0; i < spectra.size(); i++) {
138 auto &outSpec = outWS->getSpectrum(i);
139 outSpec.clearDetectorIDs();
140 for (int j = 0; j < detList[i]->xpixels(); j++)
141 for (int k = 0; k < detList[i]->ypixels(); k++)
142 outSpec.addDetectorID(static_cast<detid_t>(detList[i]->getDetectorIDAtXY(j, k)));
143 auto &outX = outSpec.mutableX();
144 auto &outY = outSpec.mutableY();
145 auto &outE = outSpec.mutableE();
146 // This is the scattered beam direction
147 V3D dir = detList[i]->getPos() - samplePos;
148
149 // Find spectra at wavelength of 1 for normalization
150 std::vector<double> xdata(1, 1.0); // wl = 1
151 std::vector<double> ydata;
152 double l2 = dir.norm();
153 // Two-theta = polar angle = scattering angle = between +Z vector and the
154 // scattered beam
155 double theta2 = dir.angle(V3D(0.0, 0.0, 1.0));
156
157 Mantid::Kernel::Unit_sptr unit = UnitFactory::Instance().create("Wavelength");
158 unit->toTOF(xdata, ydata, l1, 0,
159 {
160 {UnitParams::l2, l2},
161 {UnitParams::twoTheta, theta2},
162 });
163 double one = xdata[0];
164 double spect1 = spectrumCalc(one, iSpec, time, spectra, i);
165
166 for (size_t j = 0; j < spectra[i].size(); j++) {
167 double spect = spectra[i][j];
168
169 double relSigSpect = std::sqrt((1.0 / spect) + (1.0 / spect1));
170 if (spect1 != 0.0) {
171 spect /= spect1;
172 outX[j] = time[i][j];
173 outY[j] = spect;
174 outE[j] = relSigSpect;
175 } else {
176 throw std::runtime_error("Wavelength for normalizing to spectrum is out of range.");
177 }
178 }
179 }
180
181 Algorithm_sptr convertAlg = createChildAlgorithm("ConvertToHistogram", 0.0, 0.2);
182 convertAlg->setProperty<MatrixWorkspace_sptr>("InputWorkspace", outWS);
183 // Now execute the convert Algorithm but allow any exception to bubble up
184 convertAlg->execute();
185 outWS = convertAlg->getProperty("OutputWorkspace");
186
187 setProperty("OutputWorkspace", outWS);
188}
189
190double LoadIsawSpectrum::spectrumCalc(double TOF, int iSpec, std::vector<std::vector<double>> time,
191 std::vector<std::vector<double>> spectra, size_t id) {
192 double spect = 0;
193 if (iSpec == 1) {
194 //"Calculate the spectrum using spectral coefficients for the GSAS Type 2
195 // incident spectrum."
196 double T = TOF / 1000.; // time-of-flight in milliseconds
197
198 double c1 = spectra[id][0];
199 double c2 = spectra[id][1];
200 double c3 = spectra[id][2];
201 double c4 = spectra[id][3];
202 double c5 = spectra[id][4];
203 double c6 = spectra[id][5];
204 double c7 = spectra[id][6];
205 double c8 = spectra[id][7];
206 double c9 = spectra[id][8];
207 double c10 = spectra[id][9];
208 double c11 = spectra[id][10];
209
210 spect = c1 + c2 * exp(-c3 / std::pow(T, 2)) / std::pow(T, 5) + c4 * exp(-c5 * std::pow(T, 2)) +
211 c6 * exp(-c7 * std::pow(T, 3)) + c8 * exp(-c9 * std::pow(T, 4)) + c10 * exp(-c11 * std::pow(T, 5));
212 } else {
213 size_t i = 1;
214 for (i = 1; i < spectra[0].size() - 1; ++i)
215 if (TOF < time[id][i])
216 break;
217 spect = spectra[id][i - 1] +
218 (TOF - time[id][i - 1]) / (time[id][i] - time[id][i - 1]) * (spectra[id][i] - spectra[id][i - 1]);
219 }
220
221 return spect;
222}
223//----------------------------------------------------------------------------------------------
228 std::string grpName("Specify the Instrument");
229
230 alg->declareProperty(
231 std::make_unique<WorkspaceProperty<>>("InputWorkspace", "", Direction::Input, PropertyMode::Optional),
232 "Optional: An input workspace with the instrument we want to use.");
233
234 alg->declareProperty(std::make_unique<PropertyWithValue<std::string>>("InstrumentName", "", Direction::Input),
235 "Optional: Name of the instrument to base the "
236 "GroupingWorkspace on which to base the "
237 "GroupingWorkspace.");
238
239 alg->declareProperty(std::make_unique<FileProperty>("InstrumentFilename", "", FileProperty::OptionalLoad, ".xml"),
240 "Optional: Path to the instrument definition file on "
241 "which to base the GroupingWorkspace.");
242
243 alg->setPropertyGroup("InputWorkspace", grpName);
244 alg->setPropertyGroup("InstrumentName", grpName);
245 alg->setPropertyGroup("InstrumentFilename", grpName);
246}
247
248//----------------------------------------------------------------------------------------------
254 MatrixWorkspace_sptr inWS = alg->getProperty("InputWorkspace");
255 std::string InstrumentName = alg->getPropertyValue("InstrumentName");
256 std::string InstrumentFilename = alg->getPropertyValue("InstrumentFilename");
257
258 // Some validation
259 int numParams = 0;
260 if (inWS)
261 numParams++;
262 if (!InstrumentName.empty())
263 numParams++;
264 if (!InstrumentFilename.empty())
265 numParams++;
266
267 if (numParams > 1)
268 throw std::invalid_argument("You must specify exactly ONE way to get an "
269 "instrument (workspace, instrument name, or "
270 "IDF file). You specified more than one.");
271 if (numParams == 0)
272 throw std::invalid_argument("You must specify exactly ONE way to get an "
273 "instrument (workspace, instrument name, or "
274 "IDF file). You specified none.");
275
276 // ---------- Get the instrument one of 3 ways ---------------------------
278 if (inWS) {
279 inst = inWS->getInstrument();
280 } else {
281 Algorithm_sptr childAlg = alg->createChildAlgorithm("LoadInstrument", 0.0, 0.2);
282 MatrixWorkspace_sptr tempWS(new Workspace2D());
283 childAlg->setProperty<MatrixWorkspace_sptr>("Workspace", tempWS);
284 childAlg->setPropertyValue("Filename", InstrumentFilename);
285 childAlg->setPropertyValue("InstrumentName", InstrumentName);
286 childAlg->setProperty("RewriteSpectraMap", Mantid::Kernel::OptionalBool(false));
287 childAlg->executeAsChildAlg();
288 inst = tempWS->getInstrument();
289 }
290
291 return inst;
292}
293
294} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
Base class from which all concrete algorithm classes should be derived.
Definition: Algorithm.h:85
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
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
@ OptionalLoad
to specify a file to read but the file doesn't have to exist
Definition: FileProperty.h:53
@ Load
allowed here which will be passed to the algorithm
Definition: FileProperty.h:52
A property class for workspaces.
Load incident spectrum and detector efficiency correction file.
double spectrumCalc(double TOF, int iSpec, std::vector< std::vector< double > > time, std::vector< std::vector< double > > spectra, size_t id)
Geometry::Instrument_const_sptr getInstrument3Ways(Algorithm *alg)
Get a pointer to an instrument in one of 3 ways: InputWorkspace, InstrumentName, InstrumentFilename.
void exec() override
Run the algorithm.
void getInstrument3WaysInit(Algorithm *alg)
For use by getInstrument3Ways, initializes the properties.
Concrete workspace implementation.
Definition: Workspace2D.h:29
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
OptionalBool : Tri-state bool.
Definition: OptionalBool.h:25
The concrete, templated class for properties.
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
double angle(const V3D &) const
Angle between this and another vector.
Definition: V3D.cpp:165
double norm() const noexcept
Definition: V3D.h:263
std::shared_ptr< Algorithm > Algorithm_sptr
Typedef for a shared pointer to an Algorithm.
Definition: Algorithm.h:61
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::unique_ptr< T > create(const P &parent, const IndexArg &indexArg, const HistArg &histArg)
This is the create() method that all the other create() methods call.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::shared_ptr< Unit > Unit_sptr
Shared pointer to the Unit base class.
Definition: Unit.h:229
A namespace containing physical constants that are required by algorithms and unit routines.
Definition: Atom.h:14
int32_t detid_t
Typedef for a detector ID.
Definition: SpectrumInfo.h:21
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54