Mantid
Loading...
Searching...
No Matches
SCDPanelErrors.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 +
12#include "MantidAPI/Sample.h"
21#include <algorithm>
22#include <boost/math/special_functions/round.hpp>
23#include <cmath>
24#include <fstream>
25#include <sstream>
26#include <utility>
27
28namespace Mantid::Crystal {
29
30using namespace CurveFitting;
31
32using namespace Kernel;
33
34using namespace API;
35
36namespace {
38Logger g_log("SCDPanelErrors");
39} // namespace
40
41DECLARE_FUNCTION(SCDPanelErrors)
42
44
46SCDPanelErrors::SCDPanelErrors() : m_setupFinished(false) {
47 declareParameter("XShift", 0.0, "Shift factor in X");
48 declareParameter("YShift", 0.0, "Shift factor in Y");
49 declareParameter("ZShift", 0.0, "Shift factor in Z");
50 declareParameter("XRotate", 0.0, "Rotate angle in X");
51 declareParameter("YRotate", 0.0, "Rotate angle in Y");
52 declareParameter("ZRotate", 0.0, "Rotate angle in Z");
53 declareParameter("ScaleWidth", 1.0, "Scale width of detector");
54 declareParameter("ScaleHeight", 1.0, "Scale height of detector");
55 declareParameter("T0Shift", 0.0, "Shift for TOF");
56 declareAttribute("FileName", Attribute("", true));
57 declareAttribute("Workspace", Attribute(""));
58 declareAttribute("Bank", Attribute(""));
59}
74void SCDPanelErrors::moveDetector(double x, double y, double z, double rotx, double roty, double rotz, double scalex,
75 double scaley, std::string detname, const Workspace_sptr &inputW) const {
76 if (detname.compare("none") == 0.0)
77 return;
78 // CORELLI has sixteenpack under bank
79 DataObjects::PeaksWorkspace_sptr inputP = std::dynamic_pointer_cast<DataObjects::PeaksWorkspace>(inputW);
80 Geometry::Instrument_sptr inst = std::const_pointer_cast<Geometry::Instrument>(inputP->getInstrument());
81 if (inst->getName().compare("CORELLI") == 0.0 && detname != "moderator")
82 detname.append("/sixteenpack");
83
84 if (x != 0.0 || y != 0.0 || z != 0.0) {
85 auto alg1 = Mantid::API::AlgorithmFactory::Instance().create("MoveInstrumentComponent", -1);
86 alg1->initialize();
87 alg1->setChild(true);
88 alg1->setLogging(false);
89 alg1->setProperty<Workspace_sptr>("Workspace", inputW);
90 alg1->setPropertyValue("ComponentName", detname);
91 // Move in m
92 alg1->setProperty("X", x);
93 alg1->setProperty("Y", y);
94 alg1->setProperty("Z", z);
95 alg1->setPropertyValue("RelativePosition", "1");
96 alg1->execute();
97 }
98
99 if (rotx != 0.0) {
100 auto algx = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
101 algx->initialize();
102 algx->setChild(true);
103 algx->setLogging(false);
104 algx->setProperty<Workspace_sptr>("Workspace", inputW);
105 algx->setPropertyValue("ComponentName", detname);
106 algx->setProperty("X", 1.0);
107 algx->setProperty("Y", 0.0);
108 algx->setProperty("Z", 0.0);
109 algx->setProperty("Angle", rotx);
110 algx->setPropertyValue("RelativeRotation", "1");
111 algx->execute();
112 }
113
114 if (roty != 0.0) {
115 auto algy = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
116 algy->initialize();
117 algy->setChild(true);
118 algy->setLogging(false);
119 algy->setProperty<Workspace_sptr>("Workspace", inputW);
120 algy->setPropertyValue("ComponentName", detname);
121 algy->setProperty("X", 0.0);
122 algy->setProperty("Y", 1.0);
123 algy->setProperty("Z", 0.0);
124 algy->setProperty("Angle", roty);
125 algy->setPropertyValue("RelativeRotation", "1");
126 algy->execute();
127 }
128
129 if (rotz != 0.0) {
130 auto algz = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
131 algz->initialize();
132 algz->setChild(true);
133 algz->setLogging(false);
134 algz->setProperty<Workspace_sptr>("Workspace", inputW);
135 algz->setPropertyValue("ComponentName", detname);
136 algz->setProperty("X", 0.0);
137 algz->setProperty("Y", 0.0);
138 algz->setProperty("Z", 1.0);
139 algz->setProperty("Angle", rotz);
140 algz->setPropertyValue("RelativeRotation", "1");
141 algz->execute();
142 }
143 if (scalex != 1.0 || scaley != 1.0) {
144 Geometry::IComponent_const_sptr comp = inst->getComponentByName(detname);
145 std::shared_ptr<const Geometry::RectangularDetector> rectDet =
146 std::dynamic_pointer_cast<const Geometry::RectangularDetector>(comp);
147 if (rectDet) {
148 Geometry::ParameterMap &pmap = inputP->instrumentParameters();
149 auto oldscalex = pmap.getDouble(rectDet->getName(), "scalex");
150 auto oldscaley = pmap.getDouble(rectDet->getName(), "scaley");
151 double relscalex = scalex;
152 double relscaley = scaley;
153 if (!oldscalex.empty())
154 relscalex /= oldscalex[0];
155 if (!oldscaley.empty())
156 relscaley /= oldscaley[0];
157 pmap.addDouble(rectDet.get(), "scalex", scalex);
158 pmap.addDouble(rectDet.get(), "scaley", scaley);
159 applyRectangularDetectorScaleToComponentInfo(inputP->mutableComponentInfo(), rectDet->getComponentID(), relscalex,
160 relscaley);
161 }
162 }
163}
164
166void SCDPanelErrors::eval(double xshift, double yshift, double zshift, double xrotate, double yrotate, double zrotate,
167 double scalex, double scaley, double *out, const double *xValues, const size_t nData,
168 double tShift) const {
169 UNUSED_ARG(xValues);
170 if (nData == 0)
171 return;
172
173 setupData();
174
175 std::shared_ptr<API::Workspace> cloned = m_workspace->clone();
176 moveDetector(xshift, yshift, zshift, xrotate, yrotate, zrotate, scalex, scaley, m_bank, cloned);
177
178 auto inputP = std::dynamic_pointer_cast<DataObjects::PeaksWorkspace>(cloned);
179 // IAlgorithm_sptr alg =
180 // Mantid::API::AlgorithmFactory::Instance().create("IndexPeaks", -1);
181 // alg->initialize();
182 // alg->setChild(true);
183 // alg->setLogging(false);
184 // alg->setProperty("PeaksWorkspace", inputP);
185 // alg->setProperty("Tolerance", 0.15);
186 // alg->execute();
187 auto inst = inputP->getInstrument();
188 Geometry::OrientedLattice lattice = inputP->mutableSample().getOrientedLattice();
189 for (int i = 0; i < inputP->getNumberPeaks(); i++) {
190 const DataObjects::Peak &peak = inputP->getPeak(i);
191 V3D hkl = V3D(boost::math::iround(peak.getH()), boost::math::iround(peak.getK()), boost::math::iround(peak.getL()));
192 V3D Q2 = lattice.qFromHKL(hkl);
193 try {
194 if (hkl == V3D(0, 0, 0))
195 throw std::runtime_error("unindexed peak");
196 DataObjects::Peak peak2(inst, peak.getDetectorID(), peak.getWavelength(), hkl, peak.getGoniometerMatrix());
198
199 wl.initialize(peak2.getL1(), 0, {{UnitParams::l2, peak2.getL2()}, {UnitParams::twoTheta, peak2.getScattering()}});
200 peak2.setWavelength(wl.singleFromTOF(peak.getTOF() + tShift));
201 V3D Q3 = peak2.getQSampleFrame();
202 out[i * 3] = Q3[0] - Q2[0];
203 out[i * 3 + 1] = Q3[1] - Q2[1];
204 out[i * 3 + 2] = Q3[2] - Q2[2];
205 } catch (std::runtime_error &) {
206 // set penalty for unindexed peaks greater than tolerance
207 out[i * 3] = 0.15;
208 out[i * 3 + 1] = 0.15;
209 out[i * 3 + 2] = 0.15;
210 }
211 }
212}
213
220void SCDPanelErrors::function1D(double *out, const double *xValues, const size_t nData) const {
221 const double xshift = getParameter("XShift");
222 const double yshift = getParameter("YShift");
223 const double zshift = getParameter("ZShift");
224 const double xrotate = getParameter("XRotate");
225 const double yrotate = getParameter("YRotate");
226 const double zrotate = getParameter("ZRotate");
227 const double scalex = getParameter("ScaleWidth");
228 const double scaley = getParameter("ScaleHeight");
229 const double tShift = getParameter("T0Shift");
230 eval(xshift, yshift, zshift, xrotate, yrotate, zrotate, scalex, scaley, out, xValues, nData, tShift);
231}
232
240void SCDPanelErrors::functionDeriv1D(API::Jacobian *out, const double *xValues, const size_t nData) {
241 FunctionDomain1DView domain(xValues, nData);
242 this->calNumericalDeriv(domain, *out);
243}
244
246void SCDPanelErrors::clear() const { m_setupFinished = false; }
247
252void SCDPanelErrors::setAttribute(const std::string &attName, const IFunction::Attribute &value) {
253 if (attName == "FileName") {
254 std::string fileName = value.asUnquotedString();
255 if (fileName.empty()) {
256 storeAttributeValue("FileName", Attribute("", true));
257 return;
258 }
259 FileValidator fval;
260 std::string error = fval.isValid(fileName);
261 if (error.empty()) {
262 storeAttributeValue(attName, Attribute(fileName, true));
263 storeAttributeValue("Workspace", Attribute(""));
264 } else {
265 // file not found
266 throw Kernel::Exception::FileError(error, fileName);
267 }
268 load(fileName);
269 } else if (attName == "Workspace") {
270 std::string wsName = value.asString();
271 if (!wsName.empty()) {
272 storeAttributeValue(attName, value);
273 storeAttributeValue("FileName", Attribute("", true));
274 loadWorkspace(wsName);
275 }
276 } else {
278 m_setupFinished = false;
279 }
280}
281
286void SCDPanelErrors::load(const std::string &fname) {
287 auto loadAlg = Mantid::API::AlgorithmFactory::Instance().create("Load", -1);
288 loadAlg->initialize();
289 loadAlg->setChild(true);
290 loadAlg->setLogging(false);
291 try {
292 loadAlg->setPropertyValue("Filename", fname);
293 loadAlg->setPropertyValue("OutputWorkspace", "_SCDPanelErrors_fit_data_");
294 loadAlg->execute();
295 } catch (std::runtime_error &) {
296 throw std::runtime_error("Unable to load Nexus file for SCDPanelErrors function.");
297 }
298
299 Workspace_sptr ws = loadAlg->getProperty("OutputWorkspace");
300 API::Workspace_sptr resData = std::dynamic_pointer_cast<Mantid::API::Workspace>(ws);
301 loadWorkspace(resData);
302}
303
308void SCDPanelErrors::loadWorkspace(const std::string &wsName) const {
309 auto ws = AnalysisDataService::Instance().retrieveWS<API::Workspace>(wsName);
310 loadWorkspace(ws);
311}
312
317void SCDPanelErrors::loadWorkspace(std::shared_ptr<API::Workspace> ws) const {
318 m_workspace = std::move(ws);
319 m_setupFinished = false;
320}
321
325void SCDPanelErrors::setupData() const {
326 if (m_setupFinished) {
327 g_log.debug() << "Re-setting isn't required.";
328 return;
329 }
330
331 if (!m_workspace) {
332 std::string wsName = getAttribute("Workspace").asString();
333 if (wsName.empty())
334 throw std::invalid_argument("Data not set for function " + this->name());
335 else
336 loadWorkspace(wsName);
337 }
338
339 m_bank = getAttribute("Bank").asString();
340
341 g_log.debug() << "Setting up " << m_workspace->getName() << " bank " << m_bank << '\n';
342
343 m_setupFinished = true;
344}
345
346} // namespace Mantid::Crystal
double value
The value of the point.
Definition: FitMW.cpp:51
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
double error
Definition: IndexPeaks.cpp:133
Mantid::API::IFunction::Attribute Attribute
Definition: IsoRotDiff.cpp:25
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
1D domain - a wrapper around an array of doubles.
Attribute is a non-fitting parameter.
Definition: IFunction.h:282
void declareAttribute(const std::string &name, const API::IFunction::Attribute &defaultValue)
Declare a single attribute.
Definition: IFunction.cpp:1418
virtual void setAttribute(const std::string &name, const Attribute &)
Set a value to attribute attName.
Definition: IFunction.cpp:1409
Represents the Jacobian in IFitFunction::functionDeriv.
Definition: Jacobian.h:22
void declareParameter(const std::string &name, double initValue=0, const std::string &description="") override
Declare a new parameter.
Base Workspace Abstract Class.
Definition: Workspace.h:30
void moveDetector(double x, double y, double z, double rotx, double roty, double rotz, double scalex, double scaley, std::string detname, const API::Workspace_sptr &inputW) const
Move detectors with parameters.
std::string m_bank
Stores bank.
static const int defaultIndexValue
The default value for the workspace index.
std::shared_ptr< API::Workspace > m_workspace
Temporary workspace holder.
void setupData() const
Fill in the workspace and bank names.
void eval(double xshift, double yshift, double zshift, double xrotate, double yrotate, double zrotate, double scalex, double scaley, double *out, const double *xValues, const size_t nData, double tShift) const
Evaluate the function for a list of arguments and given scaling factor.
double getL() const override
Get the L index of the peak.
Definition: BasePeak.cpp:100
double getK() const override
Get the K index of the peak.
Definition: BasePeak.cpp:97
double getH() const override
Get the H index of the peak.
Definition: BasePeak.cpp:94
Mantid::Kernel::Matrix< double > getGoniometerMatrix() const override
Get the goniometer rotation matrix at which this peak was measured.
Definition: BasePeak.cpp:210
Structure describing a single-crystal peak.
Definition: Peak.h:34
double getL1() const override
Return the L1 flight path length (source to sample), in meters.
Definition: Peak.cpp:714
double getWavelength() const override
Calculate the neutron wavelength (in angstroms) at the peak (Note for inelastic scattering - it is th...
Definition: Peak.cpp:364
double getTOF() const override
Calculate the time of flight (in microseconds) of the neutrons for this peak, using the geometry of t...
Definition: Peak.cpp:379
int getDetectorID() const
Get the ID of the detector at the center of the peak
Definition: Peak.cpp:271
Class to implement UB matrix.
Kernel::V3D qFromHKL(const Kernel::V3D &hkl) const
Return Q-sample coordinates from hkl.
Records the filename and the description of failure.
Definition: Exception.h:98
FileValidator is a validator that checks that a filepath is valid.
Definition: FileValidator.h:25
void debug(const std::string &msg)
Logs at debug level.
Definition: Logger.cpp:114
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
void initialize(const double &_l1, const int &_emode, const UnitParametersMap &params)
Initialize the unit to perform conversion using singleToTof() and singleFromTof()
Definition: Unit.cpp:132
Wavelength in Angstrom.
Definition: Unit.h:302
Class for 3D vectors.
Definition: V3D.h:34
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
Kernel::Logger g_log("ExperimentInfo")
static logger object
MANTID_API_DLL void applyRectangularDetectorScaleToComponentInfo(Geometry::ComponentInfo &componentInfo, Geometry::IComponent *componentId, const double scaleX, const double scaleY)
Helpers for resizing RectangularDetectors.
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
Definition: IComponent.h:161
std::shared_ptr< Instrument > Instrument_sptr
Shared pointer to an instrument object.
Generate a tableworkspace to store the calibration results.
adjust instrument component position and orientation
: detector size scale at y-direction