Mantid
Loading...
Searching...
No Matches
SCDCalibratePanels2ObjFunc.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2020 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 +
7
13#include "MantidAPI/Run.h"
14#include "MantidAPI/Sample.h"
19
20#include <boost/algorithm/string.hpp>
21#include <boost/math/special_functions/round.hpp>
22#include <cmath>
23
24namespace Mantid::Crystal {
25
26using namespace Mantid::API;
27using namespace Mantid::CurveFitting;
28using namespace Mantid::DataObjects;
29using namespace Mantid::Geometry;
30using namespace Mantid::Kernel;
31
32namespace {
33// static logger
34Logger g_log("SCDCalibratePanels2ObjFunc");
35} // namespace
36
37DECLARE_FUNCTION(SCDCalibratePanels2ObjFunc)
38
39
43 // parameters for translation
44 declareParameter("DeltaX", 0.0, "relative shift along X in meter");
45 declareParameter("DeltaY", 0.0, "relative shift along Y in meter");
46 declareParameter("DeltaZ", 0.0, "relative shift along Z in meter");
47 // parameters for rotation
48 declareParameter("RotX", 0.0, "relative rotation around X in degree");
49 declareParameter("RotY", 0.0, "relative rotation around Y in degree");
50 declareParameter("RotZ", 0.0, "relative rotation around Z in degree");
51 // TOF offset for all peaks
52 // NOTE: need to have a non-zero value here
53 declareParameter("DeltaT0", 0.1, "delta of TOF");
54 // This part is for fine tuning the sample position
55 declareParameter("DeltaSampleX", 0.0, "relative shift of sample position along X.");
56 declareParameter("DeltaSampleY", 0.0, "relative shift of sample position along Y.");
57 declareParameter("DeltaSampleZ", 0.0, "relative shift of sample position along Z.");
58 // Detector size scale factors
59 declareParameter("ScaleX", 1.0, "Scale of detector along X-direction (i.e., width).");
60 declareParameter("ScaleY", 1.0, "Scale of detector along Y-direction (i.e., height).");
61}
62
63void SCDCalibratePanels2ObjFunc::setPeakWorkspace(IPeaksWorkspace_sptr &pws, const std::string &componentName,
64 const std::vector<double> &tofs) {
65 m_pws = pws->clone();
66 m_cmpt = componentName;
67
68 // Special adjustment for CORELLI
69 Instrument_sptr inst = std::const_pointer_cast<Instrument>(m_pws->getInstrument());
70 if (inst->getName().compare("CORELLI") == 0 && m_cmpt != "moderator")
71 // the second check is just to ensure that no accidental passing in
72 // a bank name with sixteenpack already appended
73 if (!boost::algorithm::ends_with(m_cmpt, "/sixteenpack"))
74 m_cmpt.append("/sixteenpack");
75
76 // Get the experimentally measured TOFs
77 m_tofs = tofs;
78
79 // Set the iteration count
80 n_iter = 0;
81}
82
91void SCDCalibratePanels2ObjFunc::function1D(double *out, const double *xValues, const size_t order) const {
92 // Get the feature vector component (numeric type)
93 //-- delta in translation
94 const double dx = getParameter("DeltaX");
95 const double dy = getParameter("DeltaY");
96 const double dz = getParameter("DeltaZ");
97 //-- delta in rotation
98 const double drx = getParameter("RotX");
99 const double dry = getParameter("RotY");
100 const double drz = getParameter("RotZ");
101 //-- delta in TOF
102 // NOTE: The T0 here is a universal offset for all peaks
103 double dT0 = getParameter("DeltaT0");
104 //-- delta of sample position
105 const double dsx = getParameter("DeltaSampleX");
106 const double dsy = getParameter("DeltaSampleY");
107 const double dsz = getParameter("DeltaSampleZ");
108 //-- scale of the detector size
109 const double scalex = getParameter("ScaleX");
110 const double scaley = getParameter("ScaleY");
111
112 //-- NOTE: given that these components are never used as
113 // one vector, there is no need to construct a
114 // xValues
115 UNUSED_ARG(xValues);
116 UNUSED_ARG(order);
117
118 // -- always working on a copy only
119 IPeaksWorkspace_sptr pws = m_pws->clone();
120
121 // NOTE: when optimizing T0, a none component will be passed in.
122 // -- For Corelli, this will be none/sixteenpack
123 // -- For others, this will be none
124 bool calibrateT0 = (m_cmpt == "none/sixteenpack") || (m_cmpt == "none");
125 // we don't need to move the instrument if we are calibrating T0
126 if (!calibrateT0) {
128
129 // translation
130 pws = moveInstruentComponentBy(dx, dy, dz, m_cmpt, pws);
131
132 // rotation
133 pws = rotateInstrumentComponentBy(drx, dry, drz, m_cmpt, pws);
134 }
135
136 // tweak sample position
137 pws = moveInstruentComponentBy(dsx, dsy, dsz, "sample-position", pws);
138
139 // calculate residual
140 // double residual = 0.0;
141 for (int i = 0; i < pws->getNumberPeaks(); ++i) {
142 // use the provided cached tofs
143 const double tof = m_tofs[i];
144
145 Peak pk = Peak(pws->getPeak(i));
146 // update instrument
147 // - this will update the instrument position attached to the peak
148 // - this will update the sample position attached to the peak
149 pk.setInstrument(pws->getInstrument());
150 // update detector ID
152 // calculate&set wavelength based on new instrument
154 wl.initialize(pk.getL1(), 0,
155 {{UnitParams::l2, pk.getL2()},
156 {UnitParams::twoTheta, pk.getScattering()},
157 {UnitParams::efixed, pk.getInitialEnergy()}});
158 pk.setWavelength(wl.singleFromTOF(tof + dT0));
159
160 V3D qv = pk.getQSampleFrame();
161 for (int j = 0; j < 3; ++j)
162 out[i * 3 + j] = qv[j];
163
164 // check the difference between n and target
165 // auto ubm = pws->sample().getOrientedLattice().getUB();
166 // V3D qv_target = ubm * pws->getPeak(i).getIntHKL();
167 // qv_target *= 2 * PI;
168 // V3D delta_qv = qv - qv_target;
169 // residual += delta_qv.norm2();
170 }
171
172 n_iter += 1;
173
174 // V3D dtrans = V3D(dx, dy, dz);
175 // V3D drots = V3D(drx, dry, drz);
176 // residual /= pws->getNumberPeaks();
177 // std::ostringstream msgiter;
178 // msgiter.precision(8);
179 // msgiter << "residual@iter_" << n_iter << ": " << residual << "\n"
180 // << "-- (dx, dy, dz) = " << dtrans << "\n"
181 // << "-- (drx, dry, drz) = " << drots << "\n"
182 // << "-- dT0 = " << dT0 << "\n\n";
183 // g_log.notice() << msgiter.str();
184}
185
186// -------///
187// Helper ///
188// -------///
189
199IPeaksWorkspace_sptr SCDCalibratePanels2ObjFunc::moveInstruentComponentBy(double deltaX, double deltaY, double deltaZ,
200 const std::string &componentName,
201 IPeaksWorkspace_sptr &pws) const {
202 // Workspace_sptr inputws = std::dynamic_pointer_cast<Workspace>(pws);
203
204 // move instrument is really fast, even with zero input
205 auto mv_alg = Mantid::API::AlgorithmFactory::Instance().create("MoveInstrumentComponent", -1);
206 mv_alg->initialize();
207 mv_alg->setChild(true);
208 mv_alg->setLogging(LOGCHILDALG);
209 mv_alg->setProperty("Workspace", pws);
210 mv_alg->setProperty("ComponentName", componentName);
211 mv_alg->setProperty("X", deltaX);
212 mv_alg->setProperty("Y", deltaY);
213 mv_alg->setProperty("Z", deltaZ);
214 mv_alg->setProperty("RelativePosition", true);
215 mv_alg->executeAsChildAlg();
216
217 return pws;
218}
219
230IPeaksWorkspace_sptr SCDCalibratePanels2ObjFunc::rotateInstrumentComponentBy(double rotX, double rotY, double rotZ,
231 const std::string &componentName,
232 IPeaksWorkspace_sptr &pws) const {
233 // rotate
234 auto rot_alg = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
235 // around X
236 rot_alg->initialize();
237 rot_alg->setChild(true);
238 rot_alg->setLogging(LOGCHILDALG);
239 rot_alg->setProperty("Workspace", pws);
240 rot_alg->setProperty("ComponentName", componentName);
241 rot_alg->setProperty("X", 1.0);
242 rot_alg->setProperty("Y", 0.0);
243 rot_alg->setProperty("Z", 0.0);
244 rot_alg->setProperty("Angle", rotX);
245 rot_alg->setProperty("RelativeRotation", true);
246 rot_alg->executeAsChildAlg();
247 // around Y
248 rot_alg->initialize();
249 rot_alg->setChild(true);
250 rot_alg->setLogging(LOGCHILDALG);
251 rot_alg->setProperty("Workspace", pws);
252 rot_alg->setProperty("ComponentName", componentName);
253 rot_alg->setProperty("X", 0.0);
254 rot_alg->setProperty("Y", 1.0);
255 rot_alg->setProperty("Z", 0.0);
256 rot_alg->setProperty("Angle", rotY);
257 rot_alg->setProperty("RelativeRotation", true);
258 rot_alg->executeAsChildAlg();
259 // around Z
260 rot_alg->initialize();
261 rot_alg->setChild(true);
262 rot_alg->setLogging(LOGCHILDALG);
263 rot_alg->setProperty("Workspace", pws);
264 rot_alg->setProperty("ComponentName", componentName);
265 rot_alg->setProperty("X", 0.0);
266 rot_alg->setProperty("Y", 0.0);
267 rot_alg->setProperty("Z", 1.0);
268 rot_alg->setProperty("Angle", rotZ);
269 rot_alg->setProperty("RelativeRotation", true);
270 rot_alg->executeAsChildAlg();
271
272 return pws;
273}
274
276SCDCalibratePanels2ObjFunc::scaleRectagularDetectorSize(const double &scalex, const double &scaley,
277 const std::string &componentName,
279
280 Geometry::Instrument_sptr inst = std::const_pointer_cast<Geometry::Instrument>(pws->getInstrument());
281 Geometry::IComponent_const_sptr comp = inst->getComponentByName(componentName);
282 std::shared_ptr<const Geometry::RectangularDetector> rectDet =
283 std::dynamic_pointer_cast<const Geometry::RectangularDetector>(comp);
284 if (rectDet) {
285 // get instrument parameter map and find out whether the
286 Geometry::ParameterMap &pmap = pws->instrumentParameters();
287 auto oldscalex = pmap.getDouble(rectDet->getName(), "scalex");
288 auto oldscaley = pmap.getDouble(rectDet->getName(), "scaley");
289 double relscalex{scalex}, relscaley{scaley};
290 if (!oldscalex.empty())
291 relscalex /= oldscalex[0];
292 if (!oldscaley.empty())
293 relscaley /= oldscaley[0];
294 applyRectangularDetectorScaleToComponentInfo(pws->mutableComponentInfo(), rectDet->getComponentID(), relscalex,
295 relscaley);
296 }
297
298 return pws;
299}
300
301} // namespace Mantid::Crystal
#define DECLARE_FUNCTION(classname)
Macro for declaring a new type of function to be used with the FunctionFactory.
#define UNUSED_ARG(x)
Function arguments are sometimes unused in certain implmentations but are required for documentation ...
Definition: System.h:64
double getParameter(size_t i) const override
Get i-th parameter.
SCDCalibratePanels2ObjFunc : TODO: DESCRIPTION.
Mantid::API::IPeaksWorkspace_sptr rotateInstrumentComponentBy(double rotX, double rotY, double rotZ, const std::string &componentName, Mantid::API::IPeaksWorkspace_sptr &pws) const
Rotate the instrument by angle axis.
void setPeakWorkspace(Mantid::API::IPeaksWorkspace_sptr &pws, const std::string &componentName, const std::vector< double > &tofs)
Mantid::API::IPeaksWorkspace_sptr scaleRectagularDetectorSize(const double &scalex, const double &scaley, const std::string &componentName, Mantid::API::IPeaksWorkspace_sptr &pws) const
void function1D(double *out, const double *xValues, const size_t order) const override
base objective function
Mantid::API::IPeaksWorkspace_sptr moveInstruentComponentBy(double deltaX, double deltaY, double deltaZ, const std::string &componentName, Mantid::API::IPeaksWorkspace_sptr &pws) const
helper functions
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
void setDetectorID(int id)
Set the detector ID of the pixel at the centre of the peak and look up and cache values related to it...
Definition: Peak.cpp:212
int getDetectorID() const
Get the ID of the detector at the center of the peak
Definition: Peak.cpp:271
void setInstrument(const Geometry::Instrument_const_sptr &inst)
Set the instrument (and save the source/sample pos).
Definition: Peak.cpp:300
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
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< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
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< 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