Mantid
Loading...
Searching...
No Matches
SCDCalibratePanels2.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
15#include "MantidAPI/Run.h"
16#include "MantidAPI/Sample.h"
17#include "MantidAPI/TableRow.h"
26#include "MantidKernel/Logger.h"
27#include <boost/container/flat_set.hpp>
28#include <boost/property_tree/ptree.hpp>
29#include <boost/property_tree/xml_parser.hpp>
30
31#include <boost/filesystem.hpp>
32#include <boost/math/special_functions/round.hpp>
33#include <cmath>
34#include <fstream>
35#include <iostream>
36#include <sstream>
37#include <utility>
38
39namespace Mantid::Crystal {
40
41using namespace Mantid::API;
42using namespace Mantid::DataObjects;
43using namespace Mantid::Geometry;
44using namespace Mantid::Kernel;
45
47namespace {
48Logger logger("SCDCalibratePanels2");
49}
50
51DECLARE_ALGORITHM(SCDCalibratePanels2)
52
53
57void SCDCalibratePanels2::init() {
58 // Input peakworkspace
59 declareProperty(std::make_unique<WorkspaceProperty<IPeaksWorkspace>>("PeakWorkspace", "", Kernel::Direction::Input),
60 "Workspace of Indexed Peaks");
61
62 // Lattice constant group
63 auto mustBeNonNegative = std::make_shared<BoundedValidator<double>>();
64 mustBeNonNegative->setLower(0.0);
65 declareProperty("RecalculateUB", true, "Recalculate UB matrix using given lattice constants");
66 declareProperty("a", EMPTY_DBL(), mustBeNonNegative,
67 "Lattice Parameter a (Leave empty to use lattice constants "
68 "in peaks workspace)");
69 declareProperty("b", EMPTY_DBL(), mustBeNonNegative,
70 "Lattice Parameter b (Leave empty to use lattice constants "
71 "in peaks workspace)");
72 declareProperty("c", EMPTY_DBL(), mustBeNonNegative,
73 "Lattice Parameter c (Leave empty to use lattice constants "
74 "in peaks workspace)");
75 declareProperty("alpha", EMPTY_DBL(), mustBeNonNegative,
76 "Lattice Parameter alpha in degrees (Leave empty to use "
77 "lattice constants in peaks workspace)");
78 declareProperty("beta", EMPTY_DBL(), mustBeNonNegative,
79 "Lattice Parameter beta in degrees (Leave empty to use "
80 "lattice constants in peaks workspace)");
81 declareProperty("gamma", EMPTY_DBL(), mustBeNonNegative,
82 "Lattice Parameter gamma in degrees (Leave empty to use "
83 "lattice constants in peaks workspace)");
84 const std::string LATTICE("Lattice Constants");
85 setPropertyGroup("RecalculateUB", LATTICE);
86 setPropertyGroup("a", LATTICE);
87 setPropertyGroup("b", LATTICE);
88 setPropertyGroup("c", LATTICE);
89 setPropertyGroup("alpha", LATTICE);
90 setPropertyGroup("beta", LATTICE);
91 setPropertyGroup("gamma", LATTICE);
92 setPropertySettings("a", std::make_unique<EnabledWhenProperty>("RecalculateUB", IS_DEFAULT));
93 setPropertySettings("b", std::make_unique<EnabledWhenProperty>("RecalculateUB", IS_DEFAULT));
94 setPropertySettings("c", std::make_unique<EnabledWhenProperty>("RecalculateUB", IS_DEFAULT));
95 setPropertySettings("alpha", std::make_unique<EnabledWhenProperty>("RecalculateUB", IS_DEFAULT));
96 setPropertySettings("beta", std::make_unique<EnabledWhenProperty>("RecalculateUB", IS_DEFAULT));
97 setPropertySettings("gamma", std::make_unique<EnabledWhenProperty>("RecalculateUB", IS_DEFAULT));
98
99 // Calibration options group
100 // NOTE:
101 // The general workflow of calibration is
102 // - calibrate L1 using all peaks
103 // - calibrate each bank
104 // - calibrate/update L1 again since bank movement will affect L1
105 // - calibrate T0
106 // - calibrate samplePos
107 const std::string CALIBRATION("Calibration Options");
108 // --------------
109 // ----- L1 -----
110 // --------------
111 declareProperty("CalibrateL1", true, "Change the L1(source to sample) distance");
112 declareProperty("SearchRadiusL1", 0.1, mustBeNonNegative,
113 "Search radius of delta L1 in meters, which is used to constrain optimization search space"
114 "when calibrating L1");
115 // editability
116 setPropertySettings("SearchRadiusL1", std::make_unique<EnabledWhenProperty>("CalibrateL1", IS_EQUAL_TO, "1"));
117 // grouping
118 setPropertyGroup("CalibrateL1", CALIBRATION);
119 setPropertyGroup("SearchRadiusL1", CALIBRATION);
120 // ----------------
121 // ----- bank -----
122 // ----------------
123 declareProperty("CalibrateBanks", false, "Calibrate position and orientation of each bank.");
124 declareProperty(
125 "SearchRadiusTransBank", 5e-2, mustBeNonNegative,
126 "This is the search radius (in meter) when calibrating component translations, used to constrain optimization"
127 "search space when calibration translation of banks");
128 declareProperty("SearchradiusRotXBank", 1.0, mustBeNonNegative,
129 "This is the search radius (in deg) when calibrating component reorientation, used to constrain "
130 "optimization search space");
131 declareProperty("SearchradiusRotYBank", 1.0, mustBeNonNegative,
132 "This is the search radius (in deg) when calibrating component reorientation, used to constrain "
133 "optimization search space");
134 declareProperty("SearchradiusRotZBank", 1.0, mustBeNonNegative,
135 "This is the search radius (in deg) when calibrating component reorientation, used to constrain "
136 "optimization search space");
137 declareProperty("CalibrateSize", false, "Calibrate detector size for each bank.");
138 declareProperty("SearchRadiusSize", 0.0, mustBeNonNegative,
139 "This is the search radius (unit less) of scale factor around at value 1.0 "
140 "when calibrating component size if it is a rectangualr detector.");
141 declareProperty("FixAspectRatio", true,
142 "If true, the scaling factor for detector along X- and Y-axis "
143 "must be the same. Otherwise, the 2 scaling factors are free.");
144 declareProperty("BankName", "",
145 "If given, only the specified bank/component will be calibrated."
146 "Otherwise, all banks will be calibrated.");
147
148 // editability
149 setPropertySettings("BankName", std::make_unique<EnabledWhenProperty>("CalibrateBanks", IS_EQUAL_TO, "1"));
150 setPropertySettings("SearchRadiusTransBank",
151 std::make_unique<EnabledWhenProperty>("CalibrateBanks", IS_EQUAL_TO, "1"));
152 setPropertySettings("SearchradiusRotXBank",
153 std::make_unique<EnabledWhenProperty>("CalibrateBanks", IS_EQUAL_TO, "1"));
154 setPropertySettings("SearchradiusRotYBank",
155 std::make_unique<EnabledWhenProperty>("CalibrateBanks", IS_EQUAL_TO, "1"));
156 setPropertySettings("SearchradiusRotZBank",
157 std::make_unique<EnabledWhenProperty>("CalibrateBanks", IS_EQUAL_TO, "1"));
158 setPropertySettings("CalibrateSize", std::make_unique<EnabledWhenProperty>("CalibrateBanks", IS_EQUAL_TO, "1"));
159 setPropertySettings("SearchRadiusSize", std::make_unique<EnabledWhenProperty>("CalibrateSize", IS_EQUAL_TO, "1"));
160 setPropertySettings("FixAspectRatio", std::make_unique<EnabledWhenProperty>("CalibrateSize", IS_EQUAL_TO, "1"));
161 // grouping
162 setPropertyGroup("CalibrateBanks", CALIBRATION);
163 setPropertyGroup("BankName", CALIBRATION);
164 setPropertyGroup("SearchRadiusTransBank", CALIBRATION);
165 setPropertyGroup("SearchradiusRotXBank", CALIBRATION);
166 setPropertyGroup("SearchradiusRotYBank", CALIBRATION);
167 setPropertyGroup("SearchradiusRotZBank", CALIBRATION);
168 setPropertyGroup("CalibrateSize", CALIBRATION);
169 setPropertyGroup("SearchRadiusSize", CALIBRATION);
170 setPropertyGroup("FixAspectRatio", CALIBRATION);
171
172 // --------------
173 // ----- T0 -----
174 // --------------
175 declareProperty("CalibrateT0", false, "Calibrate the T0 (initial TOF)");
176 declareProperty("SearchRadiusT0", 10.0, mustBeNonNegative,
177 "Search radius of T0 (in ms), used to constrain optimization search space");
178 // editability
179 setPropertySettings("SearchRadiusT0", std::make_unique<EnabledWhenProperty>("CalibrateT0", IS_EQUAL_TO, "1"));
180 // grouping
181 setPropertyGroup("CalibrateT0", CALIBRATION);
182 setPropertyGroup("SearchRadiusT0", CALIBRATION);
183 // ---------------------
184 // ----- samplePos -----
185 // ---------------------
186 declareProperty("TuneSamplePosition", false, "Fine tunning sample position");
187 declareProperty("SearchRadiusSamplePos", 0.1, mustBeNonNegative,
188 "Search radius of sample position change (in meters), used to constrain optimization search space");
189 // editability
190 setPropertySettings("SearchRadiusSamplePos",
191 std::make_unique<EnabledWhenProperty>("TuneSamplePosition", IS_EQUAL_TO, "1"));
192 // grouping
193 setPropertyGroup("TuneSamplePosition", CALIBRATION);
194 setPropertyGroup("SearchRadiusSamplePos", CALIBRATION);
195
196 // Output options group
197 declareProperty(std::make_unique<WorkspaceProperty<ITableWorkspace>>("OutputWorkspace", "", Direction::Output),
198 "The workspace containing the calibration table.");
199 declareProperty("T0", 0.0, "Returns the TOF offset from optimization", Kernel::Direction::Output);
200 const std::vector<std::string> detcalExts{".DetCal", ".Det_Cal"};
201 declareProperty(
202 std::make_unique<FileProperty>("DetCalFilename", "SCDCalibrate2.DetCal", FileProperty::OptionalSave, detcalExts),
203 "Path to an ISAW-style .detcal file to save.");
204 declareProperty(
205 std::make_unique<FileProperty>("XmlFilename", "SCDCalibrate2.xml", FileProperty::OptionalSave, ".xml"),
206 "Path to an Mantid .xml description(for LoadParameterFile) file to "
207 "save.");
208 declareProperty(
209 std::make_unique<FileProperty>("CSVFilename", "SCDCalibrate2.csv", FileProperty::OptionalSave, ".csv"),
210 "Path to an .csv file which contains the Calibration Table");
211 // group into Output group
212 const std::string OUTPUT("Output");
213 setPropertyGroup("OutputWorkspace", OUTPUT);
214 setPropertyGroup("DetCalFilename", OUTPUT);
215 setPropertyGroup("XmlFilename", OUTPUT);
216 setPropertyGroup("CSVFilename", OUTPUT);
217
218 // Add new section for advanced control of the calibration/optimization
219 // NOTE: profiling is expensive, think twice before start
220 declareProperty("VerboseOutput", false, "Toggle of child algorithm console output.");
221 declareProperty("ProfileL1", false, "Perform profiling of objective function with given input for L1");
222 declareProperty("ProfileBanks", false, "Perform profiling of objective function with given input for Banks");
223 declareProperty("ProfileT0", false, "Perform profiling of objective function with given input for T0");
224 declareProperty("ProfileL1T0", false, "Perform profiling of objective function along L1 and T0");
225 // grouping into one category
226 const std::string ADVCNTRL("Advanced Option");
227 setPropertyGroup("VerboseOutput", ADVCNTRL);
228 setPropertyGroup("ProfileL1", ADVCNTRL);
229 setPropertyGroup("ProfileBanks", ADVCNTRL);
230 setPropertyGroup("ProfileT0", ADVCNTRL);
231 setPropertyGroup("ProfileL1T0", ADVCNTRL);
232
233 // About fitting
234 auto mustBePositive = std::make_shared<Kernel::BoundedValidator<int>>();
235 mustBePositive->setLower(0);
236 declareProperty("MaxFitIterations", 500, mustBePositive,
237 "Stop after this number of iterations if a good fit is not found");
238}
239
245std::map<std::string, std::string> SCDCalibratePanels2::validateInputs() {
246 std::map<std::string, std::string> issues;
247
248 // Lattice constants are required if no UB is attached to the input
249 // peak workspace
250 IPeaksWorkspace_sptr pws = getProperty("PeakWorkspace");
251 double a = getProperty("a");
252 double b = getProperty("b");
253 double c = getProperty("c");
254 double alpha = getProperty("alpha");
255 double beta = getProperty("beta");
256 double gamma = getProperty("gamma");
257 if ((a == EMPTY_DBL() || b == EMPTY_DBL() || c == EMPTY_DBL() || alpha == EMPTY_DBL() || beta == EMPTY_DBL() ||
258 gamma == EMPTY_DBL()) &&
259 (!pws->sample().hasOrientedLattice())) {
260 issues["RecalculateUB"] = "Lattice constants are needed for peak "
261 "workspace without a UB mattrix";
262 }
263
264 // sanity check
266 throw std::runtime_error("calibrationTableColumnTypes and calibrationTableColumnTypes have different size.");
267
268 return issues;
269}
270
276 // parse all inputs
277 IPeaksWorkspace_sptr m_pws = getProperty("PeakWorkspace");
278
279 // recalculate UB with given lattice constant
280 // if required
281 if (getProperty("RecalculateUB")) {
282 // parse lattice constants
284
285 // recalculate UB and index peaks
286 updateUBMatrix(m_pws);
287 }
288
289 // remove unindexed peaks
290 m_pws = removeUnindexedPeaks(m_pws);
291
292 bool calibrateT0 = getProperty("CalibrateT0");
293 bool calibrateL1 = getProperty("CalibrateL1");
294 bool calibrateBanks = getProperty("CalibrateBanks");
295 bool tuneSamplePos = getProperty("TuneSamplePosition");
296 mCalibBankName = getPropertyValue("BankName");
297 bool profL1 = getProperty("ProfileL1");
298 bool profBanks = getProperty("ProfileBanks");
299 bool profT0 = getProperty("ProfileT0");
300 bool profL1T0 = getProperty("ProfileL1T0");
301
302 const std::string DetCalFilename = getProperty("DetCalFilename");
303 const std::string XmlFilename = getProperty("XmlFilename");
304 const std::string CSVFilename = getProperty("CSVFilename");
305
306 // Properties for resizing rectangular detector size
307 bool docalibsize = getProperty("CalibrateSize");
308 double sizesearchradius = getProperty("SearchRadiusSize");
309 bool fixdetxyratio = getProperty("FixAspectRatio");
310
311 maxFitIterations = getProperty("MaxFitIterations");
312 LOGCHILDALG = getProperty("VerboseOutput");
313
314 // STEP_0: sort the peaks
315 std::vector<std::pair<std::string, bool>> criteria{{"BankName", true}};
316 m_pws->sort(criteria);
317 // need to keep a copy of the peak workspace at its input state
318 IPeaksWorkspace_sptr pws_original = m_pws->clone();
319
320 // STEP_2: preparation
321 // get names of banks that can be calibrated
322 getBankNames(m_pws);
323
324 // DEV ONLY
325 // !!!WARNNING!!!
326 // Profiling a parameter space can be time-consuming and may freeze up your
327 // computing resources for days, therefore please proceed with caution.
328 if (profL1) {
329 profileL1(m_pws, pws_original);
330 }
331 if (profBanks) {
332 profileBanks(m_pws, pws_original);
333 }
334 if (profT0) {
335 profileT0(m_pws, pws_original);
336 }
337 if (profL1T0) {
338 profileL1T0(m_pws, pws_original);
339 }
340
341 // STEP_3: optimize
342 // - L1 (with or without T0 cali attached)
343 // - Banks
344 // - sample position
345 if (calibrateL1) {
346 // NOTE:
347 // L1 and T0 can be calibrated together to provide stable calibration results.
348 g_log.notice() << "** Calibrating L1 (moderator) as requested\n";
349 optimizeL1(m_pws, pws_original);
350 }
351
352 if (calibrateBanks) {
353 g_log.notice() << "** Calibrating L2 and orientation (bank) as requested\n";
354 optimizeBanks(m_pws, pws_original, docalibsize, sizesearchradius, fixdetxyratio);
355 }
356
357 if (calibrateL1 && calibrateBanks) {
358 g_log.notice() << "** Calibrating L1 (moderator) after bank adjusted\n";
359 optimizeL1(m_pws, pws_original);
360 // NOTE:
361 // Turns out 1 pass is sufficient (tested with the following block)
362 //
363 // double delta = 1;
364 // int cnt = 0;
365 // while (delta > 0.01) {
366 // double L1_pre = m_pws->getInstrument()->getSource()->getPos().Z();
367 // optimizeBanks(m_pws, pws_original);
368 // optimizeL1(m_pws, pws_original);
369 // double L1_post = m_pws->getInstrument()->getSource()->getPos().Z();
370 // delta = std::abs((L1_pre - L1_post) / L1_pre);
371 // cnt += 1;
372 // g_log.notice() << "@pass_" << cnt << "\n" << L1_pre << "-->" << L1_post << "\n";
373 // }
374 }
375
376 if (calibrateT0 && !calibrateL1) {
377 // NOTE:
378 // L1 and T0 can be calibrated together to provide a stable results, which is the
379 // recommended way.
380 // However, one can still calibrate T0 only if desired.
381 g_log.notice() << "** Calibrating T0 only as requested\n";
382 optimizeT0(m_pws, pws_original);
383 }
384
385 if (tuneSamplePos && !calibrateL1) {
386 g_log.notice() << "** Tunning sample position only as requested\n";
387 optimizeSamplePos(m_pws, pws_original);
388 }
389
390 if (calibrateT0 && tuneSamplePos && !calibrateL1) {
391 g_log.warning() << "** You have chosen to calibrate T0 and sample position while ignoring"
392 << " L1, which means an iterative search outside this calibration is needed"
393 << " in order to find the minimum.\n";
394 }
395
396 // STEP_4: generate a table workspace to save the calibration results
397 g_log.notice() << "-- Generate calibration table\n";
398 Instrument_sptr instCalibrated = std::const_pointer_cast<Geometry::Instrument>(m_pws->getInstrument());
399 const Geometry::ParameterMap &pmap = m_pws->instrumentParameters();
400 ITableWorkspace_sptr tablews = generateCalibrationTable(instCalibrated, pmap);
401
402 // STEP_5: Write to disk if required
403 if (!XmlFilename.empty()) {
404 saveXmlFile(XmlFilename, m_BankNames, instCalibrated, pmap);
405 }
406
407 if (!DetCalFilename.empty()) {
408 saveIsawDetCal(DetCalFilename, m_BankNames, instCalibrated, m_T0);
409 }
410
411 if (!CSVFilename.empty()) {
412 saveCalibrationTable(CSVFilename, tablews);
413 }
414
415 // STEP_4: Set the output
416 setProperty("T0", m_T0); // output the calibrated T0 as a single value
417}
418
422
430 // cache starting L1 position
431 double original_L1 = std::abs(pws->getInstrument()->getSource()->getPos().Z());
432 // T0 can be calibrate along with L1 to provide a more stable results
433 bool caliT0 = getProperty("CalibrateT0");
434 bool tuneSamplepos = getProperty("TuneSamplePosition");
435
437
438 // fit algorithm for the optimization of L1
439 auto fitL1_alg = createChildAlgorithm("Fit", -1, -1, false);
440 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
441 // NOTE: always use the original pws to get the tofs
442 std::vector<double> tofs = captureTOF(pws_original);
443 objf->setPeakWorkspace(pws, "moderator", tofs);
444 fitL1_alg->setProperty("Function", std::dynamic_pointer_cast<IFunction>(objf));
445
446 //-- bounds&constraints def
447 std::ostringstream tie_str;
448 tie_str << "DeltaX=0.0,DeltaY=0.0,"
449 << "RotX=0.0,RotY=0.0,RotZ=0.0";
450 if (!tuneSamplepos) {
451 tie_str << ",DeltaSampleX=0.0,DeltaSampleY=0.0,DeltaSampleZ=0.0";
452 }
453 if (!caliT0) {
454 tie_str << ",DeltaT0=" << m_T0;
455 }
456 std::ostringstream constraint_str;
457 double r_L1 = getProperty("SearchRadiusL1"); // get search radius
458 r_L1 = std::abs(r_L1);
459 constraint_str << -r_L1 << "<DeltaZ<" << r_L1;
460 // throw in the constrain for T0 cali if needed
461 if (caliT0) {
462 double r_dT0 = getProperty("SearchRadiusT0");
463 r_dT0 = std::abs(r_dT0);
464 constraint_str << "," << -r_dT0 << "<DeltaT0<" << r_dT0;
465 }
466 if (tuneSamplepos) {
467 double r_dsp = getProperty("SearchRadiusSamplePos");
468 r_dsp = std::abs(r_dsp);
469 constraint_str << "," << -r_dsp << "<DeltaSampleX<" << r_dsp // dsx
470 << "," << -r_dsp << "<DeltaSampleY<" << r_dsp // dsy
471 << "," << -r_dsp << "<DeltaSampleZ<" << r_dsp; // dsz
472 }
473 //-- set and go
474 fitL1_alg->setProperty("Ties", tie_str.str());
475 fitL1_alg->setProperty("Constraints", constraint_str.str());
476 fitL1_alg->setProperty("InputWorkspace", l1ws);
477 fitL1_alg->setProperty("CreateOutput", true);
478 fitL1_alg->setProperty("Output", "fit");
479 fitL1_alg->executeAsChildAlg();
480
481 //-- parse output
482 std::ostringstream calilog;
483 double chi2OverDOF = fitL1_alg->getProperty("OutputChi2overDoF");
484 ITableWorkspace_sptr rst = fitL1_alg->getProperty("OutputParameters");
485 // get results for L1
486 double dL1_optimized = rst->getRef<double>("Value", 2);
487
488 // get results for T0 (optional)
489 double dT0_optimized = rst->getRef<double>("Value", 6);
490
491 // get results for sample pos
492 // NOTE:
493 // if samplePos is not part of calibration, we will get zeros here, which means zero
494 // negative impact on the whole pws
495 double dsx_optimized = rst->getRef<double>("Value", 7);
496 double dsy_optimized = rst->getRef<double>("Value", 8);
497 double dsz_optimized = rst->getRef<double>("Value", 9);
498
499 // apply the cali results (for output cali table and file)
500 adjustComponent(0.0, 0.0, dL1_optimized, 0.0, 0.0, 0.0, EMPTY_DBL(), EMPTY_DBL(),
501 pws->getInstrument()->getSource()->getName(), pws);
502 m_T0 = dT0_optimized;
503 adjustComponent(dsx_optimized, dsy_optimized, dsz_optimized, 0.0, 0.0, 0.0, EMPTY_DBL(), EMPTY_DBL(),
504 "sample-position", pws);
505 // logging
506 int npks = pws->getNumberPeaks();
507 calilog << "-- Fit L1 results using " << npks << " peaks:\n"
508 << " dL1: " << dL1_optimized << " \n"
509 << " L1 " << original_L1 << " -> " << -pws->getInstrument()->getSource()->getPos().Z() << " \n"
510 << " dT0 = " << m_T0 << " (ms)\n"
511 << " dSamplePos = (" << dsx_optimized << "," << dsy_optimized << "," << dsz_optimized << ")\n"
512 << " chi2/DOF = " << chi2OverDOF << "\n";
513 g_log.notice() << calilog.str();
514}
515
526 const bool &docalibsize, const double &sizesearchradius,
527 const bool &fixdetxyratio) {
528
530 for (int i = 0; i < static_cast<int>(m_BankNames.size()); ++i) {
532 // prepare local copies to work with
533 const std::string bankname = *std::next(m_BankNames.begin(), i);
534 const std::string pwsBankiName = "_pws_" + bankname;
535
536 // Find out whether to skip the calibration by user's specification
537 // This check is only requied when mCalibBankName is not empty string
538 if (mCalibBankName != "") {
539 bool isbank = (bankname == mCalibBankName);
540 std::stringstream ss;
541 ss << "i = " << i << " m bank name = " << bankname;
542 if (isbank)
543 ss << " ... True ...";
544 else
545 ss << " ... Stop ...";
546 g_log.notice(ss.str());
547 // continue/skip if bank name is not what is specified
548 if (!isbank)
549 continue;
550 }
551
552 //-- step 0: extract peaks that lies on the current bank
553 IPeaksWorkspace_sptr pwsBanki = selectPeaksByBankName(pws, bankname, pwsBankiName);
554 // get tofs from the original subset of pws
555 IPeaksWorkspace_sptr pwsBanki_original = selectPeaksByBankName(pws_original, bankname, pwsBankiName);
556 std::vector<double> tofs = captureTOF(pwsBanki_original);
557
558 // Do not attempt correct panels with less than 6 peaks as the system will
559 // be under-determined
560 int nBankPeaks = pwsBanki->getNumberPeaks();
561 if (nBankPeaks < MINIMUM_PEAKS_PER_BANK) {
562 // use ostringstream to prevent OPENMP breaks log info
563 std::ostringstream msg_npeakCheckFail;
564 msg_npeakCheckFail << "-- Bank " << bankname << " have only " << nBankPeaks << " (<" << MINIMUM_PEAKS_PER_BANK
565 << ") Peaks, skipping\n";
566 g_log.notice() << msg_npeakCheckFail.str();
567 continue;
568 }
569
570 //-- step 1: prepare a mocked workspace with QSample as its yValues
572
573 //-- step 2&3: invoke fit to find both traslation and rotation
574 auto fitBank_alg = createChildAlgorithm("Fit", -1, -1, false);
575 //---- setup obj fun def
576 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
577 objf->setPeakWorkspace(pwsBanki, bankname, tofs);
578 fitBank_alg->setProperty("Function", std::dynamic_pointer_cast<IFunction>(objf));
579
580 //---- bounds&constraints def
581 //
582 double searchRadiusRotX = getProperty("SearchradiusRotXBank");
583 searchRadiusRotX = std::abs(searchRadiusRotX);
584 double searchRadiusRotY = getProperty("SearchradiusRotYBank");
585 searchRadiusRotY = std::abs(searchRadiusRotY);
586 double searchRadiusRotZ = getProperty("SearchradiusRotZBank");
587 searchRadiusRotZ = std::abs(searchRadiusRotZ);
588 //
589 double searchRadiusTran = getProperty("SearchRadiusTransBank");
590 searchRadiusTran = std::abs(searchRadiusTran);
591 //
592 std::ostringstream tie_str;
593 tie_str << "DeltaSampleX=0.0,DeltaSampleY=0.0,DeltaSampleZ=0.0,"
594 << "DeltaT0=" << m_T0;
595 std::ostringstream constraint_str;
596 // rot x
597 if (searchRadiusRotX < Tolerance) {
598 tie_str << ",RotX=0.0";
599 } else {
600 constraint_str << -searchRadiusRotX << "<RotX<" << searchRadiusRotX << ",";
601 }
602 // rot y
603 if (searchRadiusRotY < Tolerance) {
604 tie_str << ",RotY=0.0";
605 } else {
606 constraint_str << -searchRadiusRotY << "<RotY<" << searchRadiusRotY << ",";
607 }
608 // rot z
609 if (searchRadiusRotZ < Tolerance) {
610 tie_str << ",RotZ=0.0";
611 } else {
612 constraint_str << -searchRadiusRotZ << "<RotZ<" << searchRadiusRotZ << ","; // constrain rotation around Z-axis
613 }
614 // translation
615 if (searchRadiusTran < Tolerance) {
616 tie_str << ",DeltaX=0.0,DeltaY=0.0,DeltaZ=0.0";
617 } else {
618 constraint_str << -searchRadiusTran << "<DeltaX<" << searchRadiusTran << "," // restrict tranlastion along X
619 << -searchRadiusTran << "<DeltaY<" << searchRadiusTran << "," // restrict tranlastion along Y
620 << -searchRadiusTran << "<DeltaZ<" << searchRadiusTran; // restrict tranlastion along Z
621 }
622 // calibration of detector size
623 // docalibsize, sizesearchradius, fixdetxyratio
624 Geometry::Instrument_sptr inst = std::const_pointer_cast<Geometry::Instrument>(pws->getInstrument());
625 Geometry::IComponent_const_sptr comp = inst->getComponentByName(bankname);
626 std::shared_ptr<const Geometry::RectangularDetector> rectDet =
627 std::dynamic_pointer_cast<const Geometry::RectangularDetector>(comp);
628
629 std::pair<double, double> scales = getRectangularDetectorScaleFactors(inst, bankname, pws->instrumentParameters());
630
631 std::ostringstream scaleconstraints;
632 std::ostringstream scaleties;
633 if (rectDet && docalibsize) {
634 // set up constraints
635 scaleconstraints << scales.first - sizesearchradius << " <=ScaleX<" << scales.first + sizesearchradius;
636 if (fixdetxyratio) {
637 scaleties << "ScaleX=ScaleY";
638 } else {
639 scaleconstraints << "," << scales.second - sizesearchradius << " <=ScaleY<" << scales.second + sizesearchradius;
640 }
641 } else {
642 // fix the scalex and scaley to its
643 scaleties << "ScaleX=" << scales.first << ", ScaleY=" << scales.second;
644 }
645
646 // construct the final constraint and tie
647 std::string fitconstraint{constraint_str.str()};
648 if (scaleconstraints.str() != "") {
649 if (fitconstraint == "")
650 fitconstraint += scaleconstraints.str();
651 else
652 fitconstraint += "," + scaleconstraints.str();
653 }
654 std::string fittie{tie_str.str()};
655 if (scaleties.str() != "") {
656 if (fittie == "")
657 fittie += scaleties.str();
658 else
659 fittie += "," + scaleties.str();
660 }
661
662 g_log.information("Fitting " + bankname + ": constraint = " + fitconstraint + "\n\t tie = " + fittie);
663
664 //---- set&go
665 if (fittie != "")
666 fitBank_alg->setProperty("Ties", fittie);
667 if (fitconstraint != "")
668 fitBank_alg->setProperty("Constraints", fitconstraint);
669 fitBank_alg->setProperty("InputWorkspace", wsBankCali);
670 fitBank_alg->setProperty("CreateOutput", true);
671 fitBank_alg->setProperty("Output", "fit");
672 fitBank_alg->setProperty("MaxIterations", maxFitIterations);
673
674 fitBank_alg->executeAsChildAlg();
675
676 //---- cache results
677 double chi2OverDOF = fitBank_alg->getProperty("OutputChi2overDoF");
678 ITableWorkspace_sptr rstFitBank = fitBank_alg->getProperty("OutputParameters");
679 double dx = rstFitBank->getRef<double>("Value", 0);
680 double dy = rstFitBank->getRef<double>("Value", 1);
681 double dz = rstFitBank->getRef<double>("Value", 2);
682 double drx = rstFitBank->getRef<double>("Value", 3);
683 double dry = rstFitBank->getRef<double>("Value", 4);
684 double drz = rstFitBank->getRef<double>("Value", 5);
685 double scalex = rstFitBank->getRef<double>("Value", 10);
686 double scaley = rstFitBank->getRef<double>("Value", 11);
687
688 //-- step 4: update the instrument with optimization results
689 std::string bn = bankname;
690 std::ostringstream calilog;
691 if (pws->getInstrument()->getName().compare("CORELLI") == 0) {
692 bn.append("/sixteenpack");
693 }
694 // update instrument for output
695 if (rectDet && docalibsize) {
696 // adjust detector size only if it is to be set to refine
697 adjustComponent(dx, dy, dz, drx, dry, drz, scalex, scaley, bn, pws);
698 } else {
699 // (1) no rectangular det or (2) not to refine detector size:
700 // do not set any physically possible scalex or scaley
701 adjustComponent(dx, dy, dz, drx, dry, drz, EMPTY_DBL(), EMPTY_DBL(), bn, pws);
702 }
703 // logging
704 V3D dtrans(dx, dy, dz);
705 V3D drots(drx, dry, drz);
706 calilog << "-- Fit " << bn << " results using " << nBankPeaks << " peaks:\n"
707 << " d(x,y,z) = " << dtrans << "\n"
708 << " r(x,y,z) = " << drots << "\n"
709 << " scale(x, y) = " << scalex << ", " << scaley << " chi2/DOF = " << chi2OverDOF << "\n";
710 g_log.notice() << calilog.str();
711
712 // -- cleanup
714 }
716}
717
729 // create child Fit alg to optimize T0
730 auto fitT0_alg = createChildAlgorithm("Fit", -1, -1, false);
731 //-- obj func def
732 // dl;dr;
733 // Fit algorithm requires a IFunction1D to fit
734 // details
735 // Fit algorithm requires a class derived from IFunction1D as its
736 // input, so we have to implement the objective function as a separate
737 // class just to get Fit serving as an optimizer.
738 // For this particular case, we are constructing an objective function
739 // based on IFunction1D that outputs a fake histogram consist of
740 // qSample calculated based on perturbed instrument positions and
741 // orientations.
743
744 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
745 // NOTE: always use the original pws to get the tofs
746 std::vector<double> tofs = captureTOF(pws_original);
747 objf->setPeakWorkspace(pws, "none", tofs);
748 fitT0_alg->setProperty("Function", std::dynamic_pointer_cast<IFunction>(objf));
749
750 //-- bounds&constraints def
751 std::ostringstream tie_str;
752 tie_str << "DeltaX=0.0,DeltaY=0.0,DeltaZ=0.0,"
753 << "RotX=0.0,RotY=0.0,RotZ=0.0,"
754 << "DeltaSampleX=0.0,DeltaSampleY=0.0,DeltaSampleZ=0.0";
755 std::ostringstream constraint_str;
756 double r_dT0 = getProperty("SearchRadiusT0");
757 r_dT0 = std::abs(r_dT0);
758 constraint_str << -r_dT0 << "<DeltaT0<" << r_dT0;
759
760 //-- set&go
761 fitT0_alg->setProperty("Ties", tie_str.str());
762 fitT0_alg->setProperty("Constraints", constraint_str.str());
763 fitT0_alg->setProperty("InputWorkspace", t0ws);
764 fitT0_alg->setProperty("CreateOutput", true);
765 fitT0_alg->setProperty("Output", "fit");
766 fitT0_alg->executeAsChildAlg();
767
768 //-- parse output
769 std::ostringstream calilog;
770 double chi2OverDOF = fitT0_alg->getProperty("OutputChi2overDoF");
771 ITableWorkspace_sptr rst = fitT0_alg->getProperty("OutputParameters");
772 double dT0_optimized = rst->getRef<double>("Value", 6);
773
774 // apply calibration results (for output file and caliTable)
775 m_T0 = dT0_optimized;
776 int npks = pws->getNumberPeaks();
777 // logging
778 calilog << "-- Fit T0 results using " << npks << " peaks:\n"
779 << " dT0 = " << m_T0 << " (ms)\n"
780 << " chi2/DOF = " << chi2OverDOF << "\n";
781 g_log.notice() << calilog.str();
782}
783
791 // create child Fit alg to optimize T0
792 auto fitSamplePos_alg = createChildAlgorithm("Fit", -1, -1, false);
793
794 // creat input 1DHist from qSample
796
797 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
798 // NOTE: always use the original pws to get the tofs
799 std::vector<double> tofs = captureTOF(pws_original);
800 objf->setPeakWorkspace(pws, "none", tofs);
801 fitSamplePos_alg->setProperty("Function", std::dynamic_pointer_cast<IFunction>(objf));
802
803 //-- bounds&constraints def
804 std::ostringstream tie_str;
805 tie_str << "DeltaX=0.0,DeltaY=0.0,DeltaZ=0.0,"
806 << "RotX=0.0,RotY=0.0,RotZ=0.0,"
807 << "DeltaT0=" << m_T0;
808 std::ostringstream constraint_str;
809 double r_dsp = getProperty("SearchRadiusSamplePos");
810 r_dsp = std::abs(r_dsp);
811 constraint_str << -r_dsp << "<DeltaSampleX<" << r_dsp << "," << -r_dsp << "<DeltaSampleY<" << r_dsp << "," << -r_dsp
812 << "<DeltaSampleZ<" << r_dsp;
813
814 //-- set&go
815 fitSamplePos_alg->setProperty("Ties", tie_str.str());
816 fitSamplePos_alg->setProperty("Constraints", constraint_str.str());
817 fitSamplePos_alg->setProperty("InputWorkspace", samplePosws);
818 fitSamplePos_alg->setProperty("CreateOutput", true);
819 fitSamplePos_alg->setProperty("Output", "fit");
820 fitSamplePos_alg->executeAsChildAlg();
821
822 //-- parse output
823 std::ostringstream calilog;
824 double chi2OverDOF = fitSamplePos_alg->getProperty("OutputChi2overDoF");
825 ITableWorkspace_sptr rst = fitSamplePos_alg->getProperty("OutputParameters");
826 double dsx_optimized = rst->getRef<double>("Value", 7);
827 double dsy_optimized = rst->getRef<double>("Value", 8);
828 double dsz_optimized = rst->getRef<double>("Value", 9);
829
830 // apply the calibration results to pws for ouptut file
831 adjustComponent(dsx_optimized, dsy_optimized, dsz_optimized, 0.0, 0.0, 0.0, EMPTY_DBL(), EMPTY_DBL(),
832 "sample-position", pws);
833 int npks = pws->getNumberPeaks();
834 // logging
835 calilog << "-- Tune SamplePos results using " << npks << " peaks:\n"
836 << " deltaSamplePos = (" << dsx_optimized << "," << dsy_optimized << "," << dsz_optimized << ")\n"
837 << " chi2/DOF = " << chi2OverDOF << "\n";
838 g_log.notice() << calilog.str();
839}
840
844
851 m_a = getProperty("a");
852 m_b = getProperty("b");
853 m_c = getProperty("c");
854 m_alpha = getProperty("alpha");
855 m_beta = getProperty("beta");
856 m_gamma = getProperty("gamma");
857 // if any one of the six lattice constants is missing, try to get
858 // one from the workspace
859 if ((m_a == EMPTY_DBL() || m_b == EMPTY_DBL() || m_c == EMPTY_DBL() || m_alpha == EMPTY_DBL() ||
860 m_beta == EMPTY_DBL() || m_gamma == EMPTY_DBL()) &&
861 (pws->sample().hasOrientedLattice())) {
862 OrientedLattice lattice = pws->mutableSample().getOrientedLattice();
863 m_a = lattice.a();
864 m_b = lattice.b();
865 m_c = lattice.c();
866 m_alpha = lattice.alpha();
867 m_beta = lattice.beta();
868 m_gamma = lattice.gamma();
869 }
870}
871
879 auto calcUB_alg = createChildAlgorithm("CalculateUMatrix", -1, -1, false);
880 calcUB_alg->setLogging(LOGCHILDALG);
881 calcUB_alg->setProperty("PeaksWorkspace", pws);
882 calcUB_alg->setProperty("a", m_a);
883 calcUB_alg->setProperty("b", m_b);
884 calcUB_alg->setProperty("c", m_c);
885 calcUB_alg->setProperty("alpha", m_alpha);
886 calcUB_alg->setProperty("beta", m_beta);
887 calcUB_alg->setProperty("gamma", m_gamma);
888 calcUB_alg->executeAsChildAlg();
889
890 // Since UB is updated, we need to redo the indexation
891 auto idxpks_alg = createChildAlgorithm("IndexPeaks", -1, -1, false);
892 idxpks_alg->setLogging(LOGCHILDALG);
893 idxpks_alg->setProperty("PeaksWorkspace", pws);
894 idxpks_alg->setProperty("RoundHKLs", true); // both are using default
895 idxpks_alg->setProperty("Tolerance", 0.15); // values
896 idxpks_alg->executeAsChildAlg();
897}
898
906 auto fltpk_alg = createChildAlgorithm("FilterPeaks");
907 fltpk_alg->setLogging(LOGCHILDALG);
908 fltpk_alg->setProperty("InputWorkspace", pws);
909 fltpk_alg->setProperty("FilterVariable", "h^2+k^2+l^2");
910 fltpk_alg->setProperty("Operator", ">");
911 fltpk_alg->setProperty("FilterValue", 0.0);
912 fltpk_alg->setProperty("OutputWorkspace", "pws_filtered");
913 fltpk_alg->executeAsChildAlg();
914
915 IPeaksWorkspace_sptr outWS = fltpk_alg->getProperty("OutputWorkspace");
916 return outWS;
917}
918
926 std::vector<double> tofs;
927
928 for (int i = 0; i < pws->getNumberPeaks(); ++i) {
929 tofs.emplace_back(pws->getPeak(i).getTOF());
930 }
931
932 return tofs;
933}
934
941 auto peaksWorkspace = std::dynamic_pointer_cast<DataObjects::PeaksWorkspace>(pws);
942 if (!peaksWorkspace)
943 throw std::invalid_argument("a PeaksWorkspace is required to retrieve bank names");
944 int npeaks = static_cast<int>(pws->getNumberPeaks());
945 for (int i = 0; i < npeaks; ++i) {
946 std::string bname = peaksWorkspace->getPeak(i).getBankName();
947 if (bname != "None")
948 m_BankNames.insert(bname);
949 }
950}
951
961 const std::string &bankname,
962 const std::string &outputwsn) {
963 auto fltpk_alg = createChildAlgorithm("FilterPeaks");
964 fltpk_alg->setLogging(LOGCHILDALG);
965 fltpk_alg->setProperty("InputWorkspace", pws);
966 fltpk_alg->setProperty("BankName", bankname);
967 fltpk_alg->setProperty("Criterion", "=");
968 fltpk_alg->setProperty("OutputWorkspace", outputwsn);
969 fltpk_alg->executeAsChildAlg();
970
971 IPeaksWorkspace_sptr outWS = fltpk_alg->getProperty("OutputWorkspace");
972 return outWS;
973}
974
983 int npeaks = pws->getNumberPeaks();
984
985 // prepare workspace to store qSample as Histogram1D
986 MatrixWorkspace_sptr mws = std::dynamic_pointer_cast<MatrixWorkspace>(
987 WorkspaceFactory::Instance().create("Workspace2D", // use workspace 2D to mock a histogram
988 1, // one vector
989 3 * npeaks, // X :: anything is fine
990 3 * npeaks)); // Y :: flattened Q vector
991 auto &spectrum = mws->getSpectrum(0);
992 auto &xvector = spectrum.mutableX();
993 auto &yvector = spectrum.mutableY();
994 auto &evector = spectrum.mutableE();
995
996 // quick check to see what kind of weighting we can use
997 double totalSigmaInt = 0.0;
998 for (int i = 0; i < npeaks; ++i) {
999 totalSigmaInt += pws->getPeak(i).getSigmaIntensity();
1000 }
1001 double totalInt = 0.0;
1002 for (int i = 0; i < npeaks; ++i) {
1003 totalInt += pws->getPeak(i).getIntensity();
1004 }
1005 double totalCnt = 0.0;
1006 for (int i = 0; i < npeaks; ++i) {
1007 totalCnt += pws->getPeak(i).getBinCount();
1008 }
1009
1010 // directly compute qsample from UBmatrix and HKL
1011 auto ubmatrix = pws->sample().getOrientedLattice().getUB();
1012 for (int i = 0; i < npeaks; ++i) {
1013
1014 V3D qv = ubmatrix * pws->getPeak(i).getIntHKL();
1015 qv *= 2 * PI;
1016 // qv = qv / qv.norm();
1017 double wgt = 1.0;
1018 if (totalSigmaInt > 0.0) {
1019 wgt = 1.0 / pws->getPeak(i).getSigmaIntensity();
1020 } else if (totalInt > 0.0) {
1021 wgt = 1.0 / pws->getPeak(i).getIntensity();
1022 } else if (totalCnt > 0.0) {
1023 wgt = 1.0 / pws->getPeak(i).getBinCount();
1024 }
1025 // make 1dhist
1026 for (int j = 0; j < 3; ++j) {
1027 xvector[i * 3 + j] = i * 3 + j;
1028 yvector[i * 3 + j] = qv[j];
1029 evector[i * 3 + j] = wgt;
1030 }
1031 }
1032
1033 return mws;
1034}
1035
1050void SCDCalibratePanels2::adjustComponent(double dx, double dy, double dz, double drx, double dry, double drz,
1051 double scalex, double scaley, const std::string &cmptName,
1052 IPeaksWorkspace_sptr &pws) {
1053 // translation
1054 auto mv_alg = createChildAlgorithm("MoveInstrumentComponent", -1, -1, false);
1055 mv_alg->setLogging(LOGCHILDALG);
1056 mv_alg->setProperty<Workspace_sptr>("Workspace", pws);
1057 mv_alg->setProperty("ComponentName", cmptName);
1058 mv_alg->setProperty("X", dx);
1059 mv_alg->setProperty("Y", dy);
1060 mv_alg->setProperty("Z", dz);
1061 mv_alg->setProperty("RelativePosition", true);
1062 mv_alg->executeAsChildAlg();
1063
1064 // rotation
1065 auto rot_alg = createChildAlgorithm("RotateInstrumentComponent", -1, -1, false);
1066 rot_alg->setLogging(LOGCHILDALG);
1067 // - x-axis
1068 rot_alg->setProperty<Workspace_sptr>("Workspace", pws);
1069 rot_alg->setProperty("ComponentName", cmptName);
1070 rot_alg->setProperty("X", 1.0);
1071 rot_alg->setProperty("Y", 0.0);
1072 rot_alg->setProperty("Z", 0.0);
1073 rot_alg->setProperty("Angle", drx);
1074 rot_alg->setProperty("RelativeRotation", true);
1075 rot_alg->executeAsChildAlg();
1076 // - y-axis
1077 rot_alg->setProperty<Workspace_sptr>("Workspace", pws);
1078 rot_alg->setProperty("ComponentName", cmptName);
1079 rot_alg->setProperty("X", 0.0);
1080 rot_alg->setProperty("Y", 1.0);
1081 rot_alg->setProperty("Z", 0.0);
1082 rot_alg->setProperty("Angle", dry);
1083 rot_alg->setProperty("RelativeRotation", true);
1084 rot_alg->executeAsChildAlg();
1085 // - z-axis
1086 rot_alg->setProperty<Workspace_sptr>("Workspace", pws);
1087 rot_alg->setProperty("ComponentName", cmptName);
1088 rot_alg->setProperty("X", 0.0);
1089 rot_alg->setProperty("Y", 0.0);
1090 rot_alg->setProperty("Z", 1.0);
1091 rot_alg->setProperty("Angle", drz);
1092 rot_alg->setProperty("RelativeRotation", true);
1093 rot_alg->executeAsChildAlg();
1094
1095 // scale detector size
1097 auto resizeAlg = createChildAlgorithm("ResizeRectangularDetector", -1, -1, false);
1098 resizeAlg->initialize();
1099 resizeAlg->setProperty("Workspace", pws);
1100 resizeAlg->setProperty("ComponentName", cmptName);
1101 resizeAlg->setProperty("ScaleX", scalex);
1102 resizeAlg->setProperty("ScaleY", scaley);
1103 resizeAlg->execute();
1104
1105 g_log.notice() << "Resize " << cmptName << " by (absolute) " << scalex << ", " << scaley << "\n";
1106 }
1107}
1108
1116ITableWorkspace_sptr SCDCalibratePanels2::generateCalibrationTable(std::shared_ptr<Geometry::Instrument> &instrument,
1118 g_log.notice() << "Generate a TableWorkspace to store calibration results.\n";
1119
1120 // Create table workspace
1121 ITableWorkspace_sptr itablews = WorkspaceFactory::Instance().createTable();
1122 for (size_t i = 0; i < calibrationTableColumnNames.size(); ++i)
1123 itablews->addColumn(calibrationTableColumnTypes[i], calibrationTableColumnNames[i]);
1124
1125 // The first row is always the source
1126 IComponent_const_sptr source = instrument->getSource();
1127 V3D sourceRelPos = source->getRelativePos();
1128 Mantid::API::TableRow sourceRow = itablews->appendRow();
1129 // NOTE: source should not have any rotation, so we pass a zero
1130 // rotation with a fixed axis
1131 sourceRow << instrument->getSource()->getName() << sourceRelPos.X() << sourceRelPos.Y() << sourceRelPos.Z() << 1.0
1132 << 0.0 << 0.0 << 0.0 << 0.0 << 0.0;
1133
1134 // Loop through banks and set row values
1135 for (auto bankName : m_BankNames) {
1136 // CORELLLI instrument has one extra layer that pack tubes into
1137 // banks, which is what we need here
1138 if (instrument->getName().compare("CORELLI") == 0)
1139 bankName.append("/sixteenpack");
1140
1141 std::shared_ptr<const IComponent> bank = instrument->getComponentByName(bankName);
1142
1143 Quat relRot = bank->getRelativeRot();
1144 V3D pos1 = bank->getRelativePos();
1145
1146 // Calculate cosines using relRot
1147 double deg, xAxis, yAxis, zAxis;
1148 relRot.getAngleAxis(deg, xAxis, yAxis, zAxis);
1149
1150 // Detector scaling
1151 std::pair<double, double> scales = getRectangularDetectorScaleFactors(instrument, bankName, pmap);
1152
1153 // Append a new row
1154 Mantid::API::TableRow bankRow = itablews->appendRow();
1155 // Row and positions
1156 bankRow << bankName << pos1.X() << pos1.Y() << pos1.Z() << xAxis << yAxis << zAxis << deg << scales.first
1157 << scales.second;
1158 }
1159
1160 g_log.notice() << "finished generating tables\n";
1161 setProperty("OutputWorkspace", itablews);
1162
1163 return itablews;
1164}
1165
1182void SCDCalibratePanels2::saveXmlFile(const std::string &FileName,
1183 const boost::container::flat_set<std::string> &AllBankNames,
1184 std::shared_ptr<Instrument> &instrument, const Geometry::ParameterMap &pmap) {
1185 g_log.notice() << "Generating xml tree \n";
1186
1187 using boost::property_tree::ptree;
1188 ptree root;
1189 ptree parafile;
1190
1191 // configure root node
1192 parafile.put("<xmlattr>.instrument", instrument->getName());
1193 parafile.put("<xmlattr>.valid-from", instrument->getValidFromDate().toISO8601String());
1194
1195 // get L1 info for source
1196 ptree src;
1197 ptree src_dx, src_dy, src_dz;
1198 ptree src_dx_val, src_dy_val, src_dz_val;
1199 // -- get positional data from source
1200 IComponent_const_sptr source = instrument->getSource();
1201 V3D sourceRelPos = source->getRelativePos();
1202 // -- add data to node
1203 src_dx_val.put("<xmlattr>.val", sourceRelPos.X());
1204 src_dy_val.put("<xmlattr>.val", sourceRelPos.Y());
1205 src_dz_val.put("<xmlattr>.val", sourceRelPos.Z());
1206 src_dx.put("<xmlattr>.name", "x");
1207 src_dy.put("<xmlattr>.name", "y");
1208 src_dz.put("<xmlattr>.name", "z");
1209 src.put("<xmlattr>.name", source->getName());
1210
1211 src_dx.add_child("value", src_dx_val);
1212 src_dy.add_child("value", src_dy_val);
1213 src_dz.add_child("value", src_dz_val);
1214 src.add_child("parameter", src_dx);
1215 src.add_child("parameter", src_dy);
1216 src.add_child("parameter", src_dz);
1217
1218 parafile.add_child("component-link", src);
1219
1220 // add node for T0
1221 // -- property_root is a dumping group for property type values that are not tied to particular
1222 // component (i.e. virtual properties)
1223 ptree property_root;
1224 property_root.put("<xmlattr>.name", instrument->getName());
1225 ptree tof0;
1226 ptree tof0_val;
1227 tof0.put("<xmlattr>.name", "T0");
1228 tof0_val.put("<xmlattr>.val", m_T0);
1229 tof0.add_child("value", tof0_val);
1230 property_root.add_child("parameter", tof0);
1231 parafile.add_child("component-link", property_root);
1232
1233 // save sample position as a standalone component-link
1234 ptree samplePos;
1235 ptree samplePos_dx, samplePos_dy, samplePos_dz;
1236 ptree samplePos_dx_val, samplePos_dy_val, samplePos_dz_val;
1237 // -- get positional data from sample
1238 std::shared_ptr<const IComponent> sp = instrument->getComponentByName("sample-position");
1239 V3D sppos = sp->getRelativePos();
1240 samplePos_dx_val.put("<xmlattr>.val", sppos.X());
1241 samplePos_dy_val.put("<xmlattr>.val", sppos.Y());
1242 samplePos_dz_val.put("<xmlattr>.val", sppos.Z());
1243 samplePos_dx.put("<xmlattr>.name", "x");
1244 samplePos_dy.put("<xmlattr>.name", "y");
1245 samplePos_dz.put("<xmlattr>.name", "z");
1246 samplePos.put("<xmlattr>.name", "sample-position");
1247
1248 samplePos_dx.add_child("value", samplePos_dx_val);
1249 samplePos_dy.add_child("value", samplePos_dy_val);
1250 samplePos_dz.add_child("value", samplePos_dz_val);
1251 samplePos.add_child("parameter", samplePos_dx);
1252 samplePos.add_child("parameter", samplePos_dy);
1253 samplePos.add_child("parameter", samplePos_dz);
1254
1255 parafile.add_child("component-link", samplePos);
1256
1257 // configure and add each bank
1258 for (auto bankName : AllBankNames) {
1259 // Prepare data for node
1260 if (instrument->getName().compare("CORELLI") == 0)
1261 bankName.append("/sixteenpack");
1262
1263 std::shared_ptr<const IComponent> bank = instrument->getComponentByName(bankName);
1264
1265 Quat relRot = bank->getRelativeRot();
1266 std::vector<double> relRotAngles = relRot.getEulerAngles("XYZ");
1267 V3D pos1 = bank->getRelativePos();
1268 std::pair<double, double> scales = getRectangularDetectorScaleFactors(instrument, bankName, pmap);
1269
1270 // prepare node
1271 ptree bank_root;
1272 ptree bank_dx, bank_dy, bank_dz;
1273 ptree bank_dx_val, bank_dy_val, bank_dz_val;
1274 ptree bank_drotx, bank_droty, bank_drotz;
1275 ptree bank_drotx_val, bank_droty_val, bank_drotz_val;
1276 ptree bank_sx, bank_sy;
1277 ptree bank_sx_val, bank_sy_val;
1278
1279 // add data to node
1280 bank_dx_val.put("<xmlattr>.val", pos1.X());
1281 bank_dy_val.put("<xmlattr>.val", pos1.Y());
1282 bank_dz_val.put("<xmlattr>.val", pos1.Z());
1283 bank_dx.put("<xmlattr>.name", "x");
1284 bank_dy.put("<xmlattr>.name", "y");
1285 bank_dz.put("<xmlattr>.name", "z");
1286
1287 bank_drotx_val.put("<xmlattr>.val", relRotAngles[0]);
1288 bank_droty_val.put("<xmlattr>.val", relRotAngles[1]);
1289 bank_drotz_val.put("<xmlattr>.val", relRotAngles[2]);
1290 bank_drotx.put("<xmlattr>.name", "rotx");
1291 bank_droty.put("<xmlattr>.name", "roty");
1292 bank_drotz.put("<xmlattr>.name", "rotz");
1293
1294 bank_sx_val.put("<xmlattr>.val", scales.first);
1295 bank_sy_val.put("<xmlattr>.val", scales.second);
1296 bank_sx.put("<xmlattr>.name", "scalex");
1297 bank_sy.put("<xmlattr>.name", "scaley");
1298
1299 bank_root.put("<xmlattr>.name", bankName);
1300
1301 // configure structure
1302 bank_dx.add_child("value", bank_dx_val);
1303 bank_dy.add_child("value", bank_dy_val);
1304 bank_dz.add_child("value", bank_dz_val);
1305
1306 bank_drotx.add_child("value", bank_drotx_val);
1307 bank_droty.add_child("value", bank_droty_val);
1308 bank_drotz.add_child("value", bank_drotz_val);
1309
1310 bank_sx.add_child("value", bank_sx_val);
1311 bank_sy.add_child("value", bank_sy_val);
1312
1313 bank_root.add_child("parameter", bank_drotx);
1314
1315 bank_root.add_child("parameter", bank_droty);
1316 bank_root.add_child("parameter", bank_drotz);
1317 bank_root.add_child("parameter", bank_dx);
1318 bank_root.add_child("parameter", bank_dy);
1319 bank_root.add_child("parameter", bank_dz);
1320 bank_root.add_child("parameter", bank_sx);
1321 bank_root.add_child("parameter", bank_sy);
1322
1323 parafile.add_child("component-link", bank_root);
1324 }
1325
1326 // give everything to root
1327 root.add_child("parameter-file", parafile);
1328 // write the xml tree to disk
1329 g_log.notice() << "\tSaving parameter file as " << FileName << "\n";
1330 boost::property_tree::write_xml(FileName, root, std::locale(),
1331 boost::property_tree::xml_writer_settings<std::string>(' ', 2));
1332}
1333
1344void SCDCalibratePanels2::saveIsawDetCal(const std::string &filename,
1345 boost::container::flat_set<std::string> &AllBankName,
1346 std::shared_ptr<Instrument> &instrument, double T0) {
1347 g_log.notice() << "Saving DetCal file in " << filename << "\n";
1348
1349 bool tuneSamplePos = getProperty("TuneSamplePosition");
1350 if (tuneSamplePos) {
1351 g_log.warning() << "!!!WARNING!!!\n"
1352 << "DetCal format cannot retain sample position info, therefore the calibrated "
1353 << "sample position will be lost if DetCal format is the only output!\n";
1354 }
1355
1356 // create a workspace to pass to SaveIsawDetCal
1357 const size_t number_spectra = instrument->getNumberDetectors();
1358 Workspace2D_sptr wksp =
1359 std::dynamic_pointer_cast<Workspace2D>(WorkspaceFactory::Instance().create("Workspace2D", number_spectra, 2, 1));
1360 wksp->setInstrument(instrument);
1361 wksp->rebuildSpectraMapping(true /* include monitors */);
1362
1363 // convert the bank names into a vector
1364 std::vector<std::string> banknames(AllBankName.begin(), AllBankName.end());
1365
1366 // call SaveIsawDetCal
1367 auto alg = createChildAlgorithm("SaveIsawDetCal");
1368 alg->setProperty("InputWorkspace", wksp);
1369 alg->setProperty("Filename", filename);
1370 alg->setProperty("TimeOffset", T0);
1371 alg->setProperty("BankNames", banknames);
1372 alg->executeAsChildAlg();
1373}
1374
1382 auto alg = createChildAlgorithm("SaveAscii");
1383 alg->setProperty("InputWorkspace", tws);
1384 alg->setProperty("Filename", FileName);
1385 alg->setPropertyValue("CommentIndicator", "#");
1386 alg->setPropertyValue("Separator", "CSV");
1387 alg->setProperty("ColumnHeader", true);
1388 alg->setProperty("AppendToFile", false);
1389 alg->executeAsChildAlg();
1390}
1391
1399 Mantid::API::IPeaksWorkspace_sptr pws_original) {
1400 g_log.notice() << "START of profiling objective func along L1\n";
1401
1402 // control option
1403 bool verbose = getProperty("VerboseOutput");
1404 if (verbose) {
1405 // header to console
1406 g_log.notice() << "deltaL1 -- residual\n";
1407 }
1408
1409 // prepare container for profile information
1410 std::ostringstream msgrst;
1411 msgrst.precision(12);
1412 msgrst << "dL1\tresidual\n";
1413
1414 // setting up as if we are doing optimization
1415 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
1416 // NOTE: always use the original pws to get the tofs
1417 std::vector<double> tofs = captureTOF(pws_original);
1418 objf->setPeakWorkspace(pws, "moderator", tofs);
1419
1420 // call the obj to perform evaluation
1421 const int n_peaks = pws->getNumberPeaks();
1422 std::unique_ptr<double[]> target(new double[n_peaks * 3]);
1423
1424 // generate the target
1425 auto ubmatrix = pws->sample().getOrientedLattice().getUB();
1426 for (int i = 0; i < n_peaks; ++i) {
1427 V3D qv = ubmatrix * pws->getPeak(i).getIntHKL();
1428 qv *= 2 * PI;
1429 for (int j = 0; j < 3; ++j) {
1430 target[i * 3 + j] = qv[j];
1431 }
1432 }
1433
1434 double xValues[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // xValues is not used
1435
1436 // scan from -4cm to 4cm along dL1 where the minimum is supposed to be at 0 for null
1437 // case with instrument at the engineering position
1438 double deltaL1 = -4e-2;
1439 while (deltaL1 < 4e-2) {
1440 std::unique_ptr<double[]> out(new double[n_peaks * 3]);
1441 objf->setParameter("DeltaZ", deltaL1);
1442 objf->setParameter("DeltaT0", 0.0); // need to set dT0 to 0.0 if we are not cali it
1443 objf->function1D(out.get(), xValues, 1);
1444
1445 // calc residual
1446 double residual = 0.0;
1447 for (int i = 0; i < n_peaks * 3; ++i) {
1448 residual += (out[i] - target[i]) * (out[i] - target[i]);
1449 }
1450 residual = std::sqrt(residual) / (n_peaks - 1); // only 1 deg of freedom here
1451 // log rst
1452 msgrst << deltaL1 << "\t" << residual << "\n";
1453
1454 if (verbose) {
1455 g_log.notice() << deltaL1 << " -- " << residual << "\n";
1456 }
1457
1458 // increment
1459 deltaL1 += 1e-4; // 0.1mm step size
1460 }
1461
1462 // output to file
1463 auto filenamebase = boost::filesystem::temp_directory_path();
1464 filenamebase /= boost::filesystem::unique_path("profileSCDCalibratePanels2_L1.csv");
1465 std::ofstream profL1File;
1466 profL1File.open(filenamebase.string());
1467 profL1File << msgrst.str();
1468 profL1File.close();
1469 g_log.notice() << "Profile data is saved at:\n"
1470 << filenamebase << "\n"
1471 << "END of profiling objective func along L1\n";
1472}
1473
1481 const Mantid::API::IPeaksWorkspace_sptr &pws_original) {
1482 g_log.notice() << "START of profiling all banks along six degree of freedom\n";
1483
1484 // control option
1485 bool verbose = getProperty("VerboseOutput");
1486 if (verbose) {
1487 // header to console
1488 g_log.notice() << "--bankname: residual\n";
1489 }
1490
1491 // Use OPENMP to speed up the profiling
1493 for (int bankIndex = 0; bankIndex < static_cast<int>(m_BankNames.size()); ++bankIndex) {
1495 // prepare local copies to work with
1496 const std::string bankname = *std::next(m_BankNames.begin(), bankIndex);
1497 const std::string pwsBankiName = "_pws_" + bankname;
1498
1499 //-- step 0: extract peaks that lies on the current bank
1500 IPeaksWorkspace_sptr pwsBanki = selectPeaksByBankName(pws, bankname, pwsBankiName);
1501 // get tofs from the original subset of pws
1502 IPeaksWorkspace_sptr pwsBanki_original = selectPeaksByBankName(pws_original, bankname, pwsBankiName);
1503 std::vector<double> tofs = captureTOF(pwsBanki_original);
1504
1505 // Do not attempt correct panels with less than 6 peaks as the system will
1506 // be under-determined
1507 int nBankPeaks = pwsBanki->getNumberPeaks();
1508 if (nBankPeaks < MINIMUM_PEAKS_PER_BANK) {
1509 // use ostringstream to prevent OPENMP breaks log info
1510 std::ostringstream msg_npeakCheckFail;
1511 msg_npeakCheckFail << "-- Cannot profile Bank " << bankname << " have only " << nBankPeaks << " (<"
1512 << MINIMUM_PEAKS_PER_BANK << ") Peaks, skipping\n";
1513 g_log.notice() << msg_npeakCheckFail.str();
1514 continue;
1515 }
1516
1517 //
1519 std::ostringstream msgrst;
1520 msgrst.precision(12);
1521 msgrst << "dx\tdy\tdz\ttheta\tphi\trogang\tresidual\n";
1522 //
1523 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
1524 objf->setPeakWorkspace(pwsBanki, bankname, tofs);
1525 //
1526 const int n_peaks = pwsBanki->getNumberPeaks();
1527 std::unique_ptr<double[]> target(new double[n_peaks * 3]);
1528 // generate the target
1529 auto ubmatrix = pwsBanki->sample().getOrientedLattice().getUB();
1530 for (int i = 0; i < n_peaks; ++i) {
1531 V3D qv = ubmatrix * pwsBanki->getPeak(i).getIntHKL();
1532 qv *= 2 * PI;
1533 for (int j = 0; j < 3; ++j) {
1534 target[i * 3 + j] = qv[j];
1535 }
1536 }
1537 //
1538 double xValues[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // xValues is not used
1539
1540 // NOTE: very expensive scan of the parameter space
1541 for (double dx = -1e-2; dx < 1e-2; dx += 2e-2 / 20.0) {
1542 // deltaX: meter
1543 for (double dy = -1e-2; dy < 1e-2; dy += 2e-2 / 20.0) {
1544 // deltaY: meter
1545 for (double dz = -1e-2; dz < 1e-2; dz += 2e-2 / 20.0) {
1546 // deltaZ: meter
1547 for (double theta = 0.0; theta < PI; theta += PI / 20.0) {
1548 // theta: rad
1549 for (double phi = 0.0; phi < 2 * PI; phi += 2 * PI / 20.0) {
1550 // phi: rad
1551 for (double ang = -5.0; ang < 5.0; ang += 5.0 / 20.0) {
1552 // ang: degrees
1553 // configure the objfunc
1554 std::unique_ptr<double[]> out(new double[n_peaks * 3]);
1555 objf->setParameter("DeltaX", dx);
1556 objf->setParameter("DeltaY", dy);
1557 objf->setParameter("DeltaZ", dz);
1558 objf->setParameter("Theta", theta);
1559 objf->setParameter("Phi", phi);
1560 objf->setParameter("DeltaRotationAngle", ang);
1561 objf->setParameter("DeltaT0", 0.0); // need to set dT0 to 0.0 if we are not cali it
1562 objf->function1D(out.get(), xValues, 1);
1563 // calc residual
1564 double residual = 0.0;
1565 for (int i = 0; i < n_peaks * 3; ++i) {
1566 residual += (out[i] - target[i]) * (out[i] - target[i]);
1567 }
1568 residual = std::sqrt(residual) / (n_peaks - 6);
1569 // record
1570 msgrst << dx << "\t" << dy << "\t" << dz << "\t" << theta << "\t" << phi << "\t" << ang << "\t"
1571 << residual << "\n";
1572
1573 if (verbose) {
1574 g_log.notice() << "--" << bankname << ": " << residual << "\n";
1575 }
1576 }
1577 }
1578 }
1579 }
1580 }
1581 }
1582
1583 // output to file
1584 auto filenamebase = boost::filesystem::temp_directory_path();
1585 std::string fnbase = "profileSCDCalibratePanels2_" + bankname + ".csv";
1586 filenamebase /= boost::filesystem::unique_path(fnbase);
1587 std::ofstream profBankFile;
1588 profBankFile.open(filenamebase.string());
1589 profBankFile << msgrst.str();
1590 profBankFile.close();
1591
1592 // notify at the terminal
1593 std::ostringstream msg;
1594 msg << "Profile of " << bankname << " is saved at:\n"
1595 << filenamebase << "\n"
1596 << "END of profiling objective func for " << bankname << "\n";
1597 g_log.notice() << msg.str();
1599 }
1601}
1602
1610 Mantid::API::IPeaksWorkspace_sptr pws_original) {
1611 g_log.notice() << "START of profiling objective func along T0\n";
1612
1613 // control option
1614 bool verbose = getProperty("VerboseOutput");
1615 if (verbose) {
1616 // print the header to console
1617 g_log.notice() << "deltaT0 -- residual\n";
1618 }
1619
1620 // prepare container for profile information
1621 std::ostringstream msgrst;
1622 msgrst.precision(12);
1623 msgrst << "dT0\tresidual\n";
1624
1625 // setting up as if we are doing optimization
1626 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
1627 // NOTE: always use the original pws to get the tofs
1628 std::vector<double> tofs = captureTOF(pws_original);
1629 objf->setPeakWorkspace(pws, "none", tofs);
1630
1631 // generate the target
1632 const int n_peaks = pws->getNumberPeaks();
1633 std::unique_ptr<double[]> target(new double[n_peaks * 3]);
1634 auto ubmatrix = pws->sample().getOrientedLattice().getUB();
1635 for (int i = 0; i < n_peaks; ++i) {
1636 V3D qv = ubmatrix * pws->getPeak(i).getIntHKL();
1637 qv *= 2 * PI;
1638 for (int j = 0; j < 3; ++j) {
1639 target[i * 3 + j] = qv[j];
1640 }
1641 }
1642
1643 double xValues[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // xValues is not used
1644
1645 // scan from -10 ~ 10 ms along dT0
1646 double deltaT0 = -10;
1647 while (deltaT0 < 10) {
1648 std::unique_ptr<double[]> out(new double[n_peaks * 3]);
1649 objf->setParameter("DeltaT0", deltaT0);
1650 objf->function1D(out.get(), xValues, 1);
1651
1652 // calc residual
1653 double residual = 0.0;
1654 for (int i = 0; i < n_peaks * 3; ++i) {
1655 residual += (out[i] - target[i]) * (out[i] - target[i]);
1656 }
1657 residual = std::sqrt(residual) / (n_peaks - 1); // only 1 deg of freedom here
1658 // log rst
1659 msgrst << deltaT0 << "\t" << residual << "\n";
1660
1661 if (verbose) {
1662 g_log.notice() << deltaT0 << " -- " << residual << "\n";
1663 }
1664
1665 // increment
1666 deltaT0 += 0.01; // 20/2000.0
1667 }
1668
1669 // output to file
1670 auto filenamebase = boost::filesystem::temp_directory_path();
1671 filenamebase /= boost::filesystem::unique_path("profileSCDCalibratePanels2_T0.csv");
1672 std::ofstream profL1File;
1673 profL1File.open(filenamebase.string());
1674 profL1File << msgrst.str();
1675 profL1File.close();
1676 g_log.notice() << "Profile data is saved at:\n"
1677 << filenamebase << "\n"
1678 << "END of profiling objective func along T0\n";
1679}
1680
1688 Mantid::API::IPeaksWorkspace_sptr pws_original) {
1689 g_log.notice() << "START of profiling objective func along L1 and T0\n";
1690
1691 // control option
1692 bool verbose = getProperty("VerboseOutput");
1693 if (verbose) {
1694 // print the header to console
1695 g_log.notice() << "deltaL1 -- deltaT0 -- residual\n";
1696 }
1697
1698 // prepare container for profile information
1699 std::ostringstream msgrst;
1700 msgrst.precision(12);
1701 msgrst << "dL1\tdT0\tresidual\n";
1702
1703 // setting up as if we are doing optimization
1704 auto objf = std::make_shared<SCDCalibratePanels2ObjFunc>();
1705 // NOTE: always use the original pws to get the tofs
1706 std::vector<double> tofs = captureTOF(pws_original);
1707 objf->setPeakWorkspace(pws, "moderator", tofs);
1708
1709 // generate the target
1710 const int n_peaks = pws->getNumberPeaks();
1711 std::unique_ptr<double[]> target(new double[n_peaks * 3]);
1712 auto ubmatrix = pws->sample().getOrientedLattice().getUB();
1713 for (int i = 0; i < n_peaks; ++i) {
1714 V3D qv = ubmatrix * pws->getPeak(i).getIntHKL();
1715 qv *= 2 * PI;
1716 for (int j = 0; j < 3; ++j) {
1717 target[i * 3 + j] = qv[j];
1718 }
1719 }
1720
1721 double xValues[7] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; // xValues is not used
1722
1723 // profile begin
1724 for (double deltaL1 = -4e-2; deltaL1 < 4e-2; deltaL1 += 1e-4) {
1725 for (double deltaT0 = -4.0; deltaT0 < 4.0; deltaT0 += 1e-2) {
1726 std::unique_ptr<double[]> out(new double[n_peaks * 3]);
1727 objf->setParameter("DeltaZ", deltaL1);
1728 objf->setParameter("DeltaT0", deltaT0);
1729 objf->function1D(out.get(), xValues, 1);
1730
1731 // calc residual
1732 double residual = 0.0;
1733 for (int i = 0; i < n_peaks * 3; ++i) {
1734 residual += (out[i] - target[i]) * (out[i] - target[i]);
1735 }
1736 residual = std::sqrt(residual) / (n_peaks - 2); // only 1 deg of freedom here
1737
1738 if (verbose) {
1739 g_log.notice() << deltaL1 << " -- " << deltaT0 << " -- " << residual << "\n";
1740 }
1741 // log rst
1742 msgrst << deltaL1 << "\t" << deltaT0 << "\t" << residual << "\n";
1743 }
1744 }
1745
1746 // output to file
1747 auto filenamebase = boost::filesystem::temp_directory_path();
1748 filenamebase /= boost::filesystem::unique_path("profileSCDCalibratePanels2_L1T0.csv");
1749 std::ofstream profL1File;
1750 profL1File.open(filenamebase.string());
1751 profL1File << msgrst.str();
1752 profL1File.close();
1753
1754 // log
1755 g_log.notice() << "Profile data is saved at:\n"
1756 << filenamebase << "\n"
1757 << "END of profiling objective func along L1 and T0\n";
1758}
1759
1772std::pair<double, double>
1773SCDCalibratePanels2::getRectangularDetectorScaleFactors(std::shared_ptr<Geometry::Instrument> &instrument,
1774 const std::string &bankname,
1776
1777 std::pair<double, double> scales{1.0, 1.0};
1778
1779 // docalibsize, sizesearchradius, fixdetxyratio
1780 Geometry::IComponent_const_sptr comp = instrument->getComponentByName(bankname);
1781 std::shared_ptr<const Geometry::RectangularDetector> rectDet =
1782 std::dynamic_pointer_cast<const Geometry::RectangularDetector>(comp);
1783
1784 if (rectDet) {
1785 // retrieve the (scalex, scaley) stored in the workspace for this bank/component
1786 auto scalexparams = pmap.getDouble(rectDet->getName(), "scalex");
1787 auto scaleyparams = pmap.getDouble(rectDet->getName(), "scaley");
1788 if (!scalexparams.empty())
1789 scales.first = scalexparams[0];
1790 if (!scaleyparams.empty())
1791 scales.second = scaleyparams[0];
1792 }
1793
1794 return scales;
1795}
1796
1797} // namespace Mantid::Crystal
#define DECLARE_ALGORITHM(classname)
Definition: Algorithm.h:576
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
Definition: MultiThreaded.h:94
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
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
Kernel::Logger & g_log
Definition: Algorithm.h:451
@ OptionalSave
to specify a file to write to but an empty string is
Definition: FileProperty.h:50
TableRow represents a row in a TableWorkspace.
Definition: TableRow.h:39
A property class for workspaces.
SCDCalibratePanels2 : Using input peakworkspace with indexation results to calibrate each individual ...
Mantid::API::IPeaksWorkspace_sptr removeUnindexedPeaks(const Mantid::API::IPeaksWorkspace_sptr &pws)
Remove unindexed peaks from workspace.
void adjustComponent(double dx, double dy, double dz, double drx, double dry, double drz, double scalex, double scaley, const std::string &cmptName, Mantid::API::IPeaksWorkspace_sptr &pws)
Helper functions for adjusting components.
void saveIsawDetCal(const std::string &filename, boost::container::flat_set< std::string > &AllBankName, std::shared_ptr< Geometry::Instrument > &instrument, double T0)
Save to ISAW type det calibration output for backward compatiblity.
Mantid::API::ITableWorkspace_sptr generateCalibrationTable(std::shared_ptr< Geometry::Instrument > &instrument, const Geometry::ParameterMap &pmap)
Generate a Table workspace to store the calibration results.
void optimizeT0(Mantid::API::IPeaksWorkspace_sptr pws, Mantid::API::IPeaksWorkspace_sptr pws_original)
Private function for calibrating T0.
void optimizeL1(Mantid::API::IPeaksWorkspace_sptr pws, Mantid::API::IPeaksWorkspace_sptr pws_original)
Private function for calibrating L1.
double m_a
unique vars for a given instance of calibration
void optimizeSamplePos(Mantid::API::IPeaksWorkspace_sptr pws, Mantid::API::IPeaksWorkspace_sptr pws_original)
Private function for fine tunning sample position.
const std::vector< std::string > calibrationTableColumnTypes
void getBankNames(const Mantid::API::IPeaksWorkspace_sptr &pws)
Private function for getting names of banks to be calibrated.
void profileBanks(Mantid::API::IPeaksWorkspace_sptr &pws, const Mantid::API::IPeaksWorkspace_sptr &pws_original)
Profiling obj func along six degree of freedom, which can be very slow.
void saveCalibrationTable(const std::string &FileName, Mantid::API::ITableWorkspace_sptr &tws)
Save the calibration table to a CSV file.
boost::container::flat_set< std::string > m_BankNames
void updateUBMatrix(const Mantid::API::IPeaksWorkspace_sptr &pws)
Update the UB matrix.
void optimizeBanks(Mantid::API::IPeaksWorkspace_sptr pws, const Mantid::API::IPeaksWorkspace_sptr &pws_original, const bool &docalibsize, const double &sizesearchradius, const bool &fixdetxyratio)
Private function for calibrating banks.
void parseLatticeConstant(const Mantid::API::IPeaksWorkspace_sptr &pws)
Private function dedicated for parsing lattice constant.
const std::vector< std::string > calibrationTableColumnNames
void saveXmlFile(const std::string &FileName, const boost::container::flat_set< std::string > &AllBankNames, std::shared_ptr< Geometry::Instrument > &instrument, const Geometry::ParameterMap &pmap)
Save to xml file for Mantid to load by manual crafting.
void exec() override
Overwrites Algorithm method.
void profileT0(Mantid::API::IPeaksWorkspace_sptr &pws, Mantid::API::IPeaksWorkspace_sptr pws_original)
Profile obj func along T0 axis.
Mantid::API::IPeaksWorkspace_sptr selectPeaksByBankName(const Mantid::API::IPeaksWorkspace_sptr &pws, const std::string &bankname, const std::string &outputwsn)
Helper function for selecting peaks based on given bank name.
std::map< std::string, std::string > validateInputs() override
Private validator for inputs.
std::pair< double, double > getRectangularDetectorScaleFactors(std::shared_ptr< Geometry::Instrument > &instrument, const std::string &bankname, const Geometry::ParameterMap &pmap)
Retrieve "scalex" and "scaley" from a workspace's parameter map if the component is rectangular detec...
Mantid::API::MatrixWorkspace_sptr getIdealQSampleAsHistogram1D(const Mantid::API::IPeaksWorkspace_sptr &pws)
Helper function that calculates the ideal qSample based on integer HKL.
void profileL1(Mantid::API::IPeaksWorkspace_sptr &pws, Mantid::API::IPeaksWorkspace_sptr pws_original)
Profile related functions.
void profileL1T0(Mantid::API::IPeaksWorkspace_sptr &pws, Mantid::API::IPeaksWorkspace_sptr pws_original)
Profile obj func along L1 and T0 axis.
std::vector< double > captureTOF(const Mantid::API::IPeaksWorkspace_sptr &pws)
Cache TOF equivalent to those measured from experiment.
Class to implement UB matrix.
double alpha() const
Get lattice parameter.
Definition: UnitCell.cpp:133
double a(int nd) const
Get lattice parameter a1-a3 as function of index (0-2)
Definition: UnitCell.cpp:94
double c() const
Get lattice parameter.
Definition: UnitCell.cpp:128
double beta() const
Get lattice parameter.
Definition: UnitCell.cpp:138
double b() const
Get lattice parameter.
Definition: UnitCell.cpp:123
double gamma() const
Get lattice parameter.
Definition: UnitCell.cpp:143
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
The Logger class is in charge of the publishing messages from the framework through various channels.
Definition: Logger.h:52
void notice(const std::string &msg)
Logs at notice level.
Definition: Logger.cpp:95
void warning(const std::string &msg)
Logs at warning level.
Definition: Logger.cpp:86
void information(const std::string &msg)
Logs at information level.
Definition: Logger.cpp:105
Class for quaternions.
Definition: Quat.h:39
void getAngleAxis(double &_deg, double &_ax0, double &_ax1, double &ax2) const
Extracts the angle of roatation and axis.
Definition: Quat.cpp:135
std::vector< double > getEulerAngles(const std::string &convention) const
Calculate the Euler angles that are equivalent to this Quaternion.
Definition: Quat.cpp:729
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
constexpr double X() const noexcept
Get x.
Definition: V3D.h:232
constexpr double Y() const noexcept
Get y.
Definition: V3D.h:233
constexpr double Z() const noexcept
Get z.
Definition: V3D.h:234
std::shared_ptr< IPeaksWorkspace > IPeaksWorkspace_sptr
shared pointer to Mantid::API::IPeaksWorkspace
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< Workspace > Workspace_sptr
shared pointer to Mantid::API::Workspace
Definition: Workspace_fwd.h:20
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
std::shared_ptr< Workspace2D > Workspace2D_sptr
shared pointer to Mantid::DataObjects::Workspace2D
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 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.
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
Definition: MultiThreaded.h:22
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition: EmptyValues.h:43
Generate a tableworkspace to store the calibration results.
adjust instrument component position and orientation
: detector size scale at y-direction
@ Input
An input workspace.
Definition: Property.h:53
@ Output
An output workspace.
Definition: Property.h:54