39using namespace Kernel;
41using namespace Geometry;
42using namespace DataObjects;
43using Types::Core::DateAndTime;
52 double x,
y,
z, rotx, roty, rotz;
53 std::string detname, inname, outname, peakOpt, rb_param, groupWSName;
54 auto *p =
reinterpret_cast<std::string *
>(params);
61 x = gsl_vector_get(v, 0);
62 y = gsl_vector_get(v, 1);
63 z = gsl_vector_get(v, 2);
64 rotx = gsl_vector_get(v, 3);
65 roty = gsl_vector_get(v, 4);
66 rotz = gsl_vector_get(v, 5);
68 return u.
intensity(
x,
y,
z, rotx, roty, rotz, detname, inname, outname, peakOpt, rb_param, groupWSName);
84 double rotz,
const std::string &detname,
89 alg1->setPropertyValue(
"ComponentName", detname);
91 alg1->setProperty(
"X",
x * 0.01);
92 alg1->setProperty(
"Y",
y * 0.01);
93 alg1->setProperty(
"Z",
z * 0.01);
94 alg1->setPropertyValue(
"RelativePosition",
"1");
95 alg1->executeAsChildAlg();
99 algx->setPropertyValue(
"ComponentName", detname);
100 algx->setProperty(
"X", 1.0);
101 algx->setProperty(
"Y", 0.0);
102 algx->setProperty(
"Z", 0.0);
103 algx->setProperty(
"Angle", rotx);
104 algx->setPropertyValue(
"RelativeRotation",
"1");
105 algx->executeAsChildAlg();
109 algy->setPropertyValue(
"ComponentName", detname);
110 algy->setProperty(
"X", 0.0);
111 algy->setProperty(
"Y", 1.0);
112 algy->setProperty(
"Z", 0.0);
113 algy->setProperty(
"Angle", roty);
114 algy->setPropertyValue(
"RelativeRotation",
"1");
115 algy->executeAsChildAlg();
119 algz->setPropertyValue(
"ComponentName", detname);
120 algz->setProperty(
"X", 0.0);
121 algz->setProperty(
"Y", 0.0);
122 algz->setProperty(
"Z", 1.0);
123 algz->setProperty(
"Angle", rotz);
124 algz->setPropertyValue(
"RelativeRotation",
"1");
125 algz->executeAsChildAlg();
145 double rotz,
const std::string &detname,
const std::string &inname,
146 const std::string &outname,
const std::string &peakOpt,
147 const std::string &rb_param,
const std::string &groupWSName) {
155 g_log.
debug() << tim <<
" to movedetector()\n";
159 alg3->setPropertyValue(
"OutputWorkspace", outname);
160 alg3->setPropertyValue(
"Target",
"dSpacing");
161 alg3->executeAsChildAlg();
169 alg4->setPropertyValue(
"GroupingFileName",
"");
170 alg4->setPropertyValue(
"GroupingWorkspace", groupWSName);
171 alg4->executeAsChildAlg();
172 outputW = alg4->getProperty(
"OutputWorkspace");
175 g_log.
debug() << tim <<
" to DiffractionFocussing\n";
180 alg5->setPropertyValue(
"Params", rb_param);
181 alg5->executeAsChildAlg();
182 outputW = alg5->getProperty(
"OutputWorkspace");
187 const MantidVec &yValues = outputW->readY(0);
188 auto it = std::max_element(yValues.begin(), yValues.end());
189 double peakHeight = *it;
192 double peakLoc = outputW->readX(0)[it - yValues.begin()];
202 std::ostringstream fun_str;
203 fun_str <<
"name=Gaussian,Height=" << peakHeight <<
",Sigma=0.01,PeakCentre=" << peakLoc;
204 fit_alg->setProperty(
"Function", fun_str.str());
205 fit_alg->setProperty(
"InputWorkspace", outputW);
206 fit_alg->setProperty(
"WorkspaceIndex", 0);
207 fit_alg->setProperty(
"StartX", outputW->readX(0)[0]);
208 fit_alg->setProperty(
"EndX", outputW->readX(0)[outputW->blocksize()]);
209 fit_alg->setProperty(
"MaxIterations", 200);
210 fit_alg->setProperty(
"Output",
"fit");
211 fit_alg->executeAsChildAlg();
215 std::vector<double> params;
217 for (
size_t i = 0; i < fun_res->nParams(); ++i) {
218 params.emplace_back(fun_res->getParameter(i));
220 peakHeight = params[0];
225 g_log.
debug() << tim <<
" to movedetector()\n";
230 return (
static_cast<int>(inputE->getNumberEvents()) / 1.e6) / peakHeight +
231 std::fabs(peakLoc - boost::lexical_cast<double>(peakOpt));
238 std::make_shared<InstrumentValidator>()),
239 "The workspace containing the geometry to be calibrated.");
242 "A comma separated list of first bin boundary, width, last "
243 "bin boundary. Optionally "
244 "this can be followed by a comma and more widths and last "
246 "Use bin boundaries close to peak you wish to maximize. "
247 "Negative width values indicate logarithmic binning.");
249 auto mustBePositive = std::make_shared<BoundedValidator<int>>();
251 "Stop after this number of iterations if a good fit is not found");
253 auto dblmustBePositive = std::make_shared<BoundedValidator<double>>();
255 "Optimize this location of peak by moving detectors");
258 "The output filename of the ISAW DetCal file");
261 "Optional: To only calibrate one bank. Any bank whose name does not "
262 "match the given string will have no events.");
265 gsl_set_error_handler_off();
274 const int maxIterations =
getProperty(
"MaxIterations");
275 const double peakOpt =
getProperty(
"LocationOfPeakToOptimize");
281 const std::string rb_params =
getProperty(
"Params");
286 const auto &dummyW = create<EventWorkspace>(*inputW, 1, inputW->binEdges(0));
290 std::vector<std::shared_ptr<RectangularDetector>> detList;
293 bool doOneBank = (!onebank.empty());
294 for (
int i = 0; i < inst->nelements(); i++) {
295 std::shared_ptr<RectangularDetector> det;
296 std::shared_ptr<ICompAssembly> assem;
297 std::shared_ptr<ICompAssembly> assem2;
299 det = std::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
301 if (det->getName() == onebank)
302 detList.emplace_back(det);
304 detList.emplace_back(det);
309 assem = std::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
311 for (
int j = 0; j < assem->nelements(); j++) {
312 det = std::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
314 if (det->getName() == onebank)
315 detList.emplace_back(det);
317 detList.emplace_back(det);
324 assem2 = std::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
326 for (
int k = 0; k < assem2->nelements(); k++) {
327 det = std::dynamic_pointer_cast<RectangularDetector>((*assem2)[k]);
329 if (det->getName() == onebank)
330 detList.emplace_back(det);
332 detList.emplace_back(det);
344 std::string inname =
getProperty(
"InputWorkspace");
345 std::string outname = inname +
"2";
348 algS->setProperty(
"InputWorkspace", inputW);
349 algS->setPropertyValue(
"SortBy",
"X Value");
350 algS->executeAsChildAlg();
353 std::string filename =
getProperty(
"DetCalFilename");
354 std::fstream outfile;
355 outfile.open(filename.c_str(), std::ios::out);
357 if (detList.size() > 1) {
359 outfile <<
"# Mantid Optimized .DetCal file for SNAP with TWO detector "
361 outfile <<
"# Old Panel, nominal size and distance at -90 degrees.\n";
362 outfile <<
"# New Panel, nominal size and distance at +90 degrees.\n";
364 outfile <<
"# Lengths are in centimeters.\n";
365 outfile <<
"# Base and up give directions of unit vectors for a local\n";
366 outfile <<
"# x,y coordinate system on the face of the detector.\n";
368 outfile <<
"# " << DateAndTime::getCurrentTime().toFormattedString(
"%c") <<
"\n";
370 outfile <<
"6 L1 T0_SHIFT\n";
373 outfile <<
"7 " << source->getDistance(*sample) * 100 <<
" 0\n";
374 outfile <<
"4 DETNUM NROWS NCOLS WIDTH HEIGHT DEPTH DETD "
375 "CenterX CenterY CenterZ BaseX BaseY BaseZ "
379 Progress prog(
this, 0.0, 1.0, detList.size());
380 for (
int det = 0; det < static_cast<int>(detList.size()); det++) {
382 par[0] = detList[det]->getName();
385 std::ostringstream strpeakOpt;
386 strpeakOpt << peakOpt;
387 par[3] = strpeakOpt.str();
394 alg2->setProperty(
"InputWorkspace", inputW);
395 alg2->setPropertyValue(
"GroupNames", detList[det]->
getName());
396 std::string groupWSName =
"group_" + detList[det]->getName();
397 alg2->setPropertyValue(
"OutputWorkspace", groupWSName);
398 alg2->executeAsChildAlg();
399 par[5] = groupWSName;
400 std::cout << tim <<
" to CreateGroupingWorkspace\n";
402 const gsl_multimin_fminimizer_type *T = gsl_multimin_fminimizer_nmsimplex;
404 gsl_multimin_function minex_func;
413 x = gsl_vector_alloc(nopt);
414 gsl_vector_set(
x, 0, 0.0);
415 gsl_vector_set(
x, 1, 0.0);
416 gsl_vector_set(
x, 2, 0.0);
417 gsl_vector_set(
x, 3, 0.0);
418 gsl_vector_set(
x, 4, 0.0);
419 gsl_vector_set(
x, 5, 0.0);
422 ss = gsl_vector_alloc(nopt);
423 gsl_vector_set_all(ss, 0.1);
428 minex_func.params = ∥
430 gsl_multimin_fminimizer *s = gsl_multimin_fminimizer_alloc(T, nopt);
431 gsl_multimin_fminimizer_set(s, &minex_func,
x, ss);
435 status = gsl_multimin_fminimizer_iterate(s);
440 double size = gsl_multimin_fminimizer_size(s);
441 status = gsl_multimin_test_size(size, 1e-2);
443 }
while (status == GSL_CONTINUE && iter < maxIterations && s->fval != -0.000);
446 if (s->fval != -0.000)
447 movedetector(gsl_vector_get(s->x, 0), gsl_vector_get(s->x, 1), gsl_vector_get(s->x, 2), gsl_vector_get(s->x, 3),
448 gsl_vector_get(s->x, 4), gsl_vector_get(s->x, 5), par[0],
getProperty(
"InputWorkspace"));
450 gsl_vector_set(s->x, 0, 0.0);
451 gsl_vector_set(s->x, 1, 0.0);
452 gsl_vector_set(s->x, 2, 0.0);
453 gsl_vector_set(s->x, 3, 0.0);
454 gsl_vector_set(s->x, 4, 0.0);
455 gsl_vector_set(s->x, 5, 0.0);
458 std::string reportOfDiffractionEventCalibrateDetectors = gsl_strerror(status);
459 if (s->fval == -0.000)
460 reportOfDiffractionEventCalibrateDetectors =
"No events";
466 <<
"Iteration = " << iter <<
"\n"
467 <<
"Status = " << reportOfDiffractionEventCalibrateDetectors <<
"\n"
468 <<
"Minimize PeakLoc-" << peakOpt <<
" = " << s->fval <<
"\n";
470 g_log.
information() <<
"Move (X) = " << gsl_vector_get(s->x, 0) * 0.01 <<
" \n";
471 g_log.
information() <<
"Move (Y) = " << gsl_vector_get(s->x, 1) * 0.01 <<
" \n";
472 g_log.
information() <<
"Move (Z) = " << gsl_vector_get(s->x, 2) * 0.01 <<
" \n";
478 V3D(gsl_vector_get(s->x, 0) * 0.01, gsl_vector_get(s->x, 1) * 0.01, gsl_vector_get(s->x, 2) * 0.01);
479 Kernel::V3D Center = detList[det]->getPos() + CalCenter;
480 int pixmax = detList[det]->xpixels() - 1;
481 int pixmid = (detList[det]->ypixels() - 1) / 2;
483 detList[det]->getAtXY(pixmax, pixmid)->getBoundingBox(box);
484 double baseX = box.
xMax();
485 double baseY = box.
yMax();
486 double baseZ = box.
zMax();
488 pixmid = (detList[det]->xpixels() - 1) / 2;
489 pixmax = detList[det]->ypixels() - 1;
490 detList[det]->getAtXY(pixmid, pixmax)->getBoundingBox(box);
491 double upX = box.
xMax();
492 double upY = box.
yMax();
493 double upZ = box.
zMax();
502 double angle = gsl_vector_get(s->x, 3) *
deg2rad;
503 Base =
V3D(baseX, baseY * cos(angle) - baseZ * sin(angle), baseY * sin(angle) + baseZ * cos(angle));
507 Up =
V3D(upX, upY * cos(angle) - upZ * sin(angle), upY * sin(angle) + upZ * cos(angle));
512 angle = gsl_vector_get(s->x, 4) *
deg2rad;
513 Base =
V3D(baseZ * sin(angle) + baseX * cos(angle), baseY, baseZ * cos(angle) - baseX * sin(angle));
517 Up =
V3D(upZ * cos(angle) - upX * sin(angle), upY, upZ * sin(angle) + upX * cos(angle));
522 angle = gsl_vector_get(s->x, 5) *
deg2rad;
523 Base =
V3D(baseX * cos(angle) - baseY * sin(angle), baseX * sin(angle) + baseY * cos(angle), baseZ);
527 Up =
V3D(upX * cos(angle) - upY * sin(angle), upX * sin(angle) + upY * cos(angle), upZ);
532 outfile <<
"5 " << detList[det]->getName().substr(4) <<
" " << detList[det]->xpixels() <<
" "
533 << detList[det]->ypixels() <<
" " << 100.0 * detList[det]->xsize() <<
" " << 100.0 * detList[det]->ysize()
536 <<
" " << Center.
norm() <<
" ";
537 Center.
write(outfile);
547 gsl_multimin_fminimizer_free(s);
#define DECLARE_ALGORITHM(classname)
std::string getName(const IMDDimension &self)
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
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.
@ Save
to specify a file to write to, the file may or may not exist
Helper class for reporting progress from algorithms.
A property class for workspaces.
Find the offsets for each detector.
void movedetector(double x, double y, double z, double rotx, double roty, double rotz, const std::string &detname, const Mantid::DataObjects::EventWorkspace_sptr &inputW)
The movedetector function changes detector position and angles.
void init() override
Initialisation method.
double intensity(double x, double y, double z, double rotx, double roty, double rotz, const std::string &detname, const std::string &inname, const std::string &outname, const std::string &peakOpt, const std::string &rb_param, const std::string &groupWSName)
Function to optimize.
void exec() override
Executes the algorithm.
A simple structure that defines an axis-aligned cuboid shaped bounding box for a geometrical object.
double xMax() const
Return the maximum value of X.
double zMax() const
Return the maximum value of Z.
double yMax() const
Return the maximum value of Y.
CPUTimer : Timer that uses the CPU time, rather than wall-clock time to measure execution time.
Exception for when an item is not found in a collection.
void debug(const std::string &msg)
Logs at debug level.
void error(const std::string &msg)
Logs at error 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.
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...
double normalize()
Make a normalized vector (return norm value)
double norm() const noexcept
void write(std::ostream &) const
Write out the point values.
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< IFunction > IFunction_sptr
shared pointer to the function base class
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
static double gsl_costFunction(const gsl_vector *v, void *params)
The gsl_costFunction is optimized by GSL simplex.
std::shared_ptr< const EventWorkspace > EventWorkspace_const_sptr
shared pointer to a const Workspace2D
std::shared_ptr< EventWorkspace > EventWorkspace_sptr
shared pointer to the EventWorkspace class
constexpr double deg2rad
Defines units/enum for Crystal work.
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::vector< double > MantidVec
typedef for the data storage used in Mantid matrix workspaces
@ Input
An input workspace.