Mantid
Loading...
Searching...
No Matches
SCDCalibratePanels.cpp
Go to the documentation of this file.
1// Mantid Repository : https://github.com/mantidproject/mantid
2//
3// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
4// NScD Oak Ridge National Laboratory, European Spallation Source,
5// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
6// SPDX - License - Identifier: GPL - 3.0 +
13#include "MantidAPI/IFunction.h"
15#include "MantidAPI/Run.h"
16#include "MantidAPI/Sample.h"
30#include <boost/container/flat_set.hpp>
31#include <boost/math/special_functions/round.hpp>
32#include <fstream>
33#include <sstream>
34
35using namespace Mantid::DataObjects;
36using namespace Mantid::API;
37using namespace std;
38using namespace Mantid::Geometry;
39using namespace Mantid::Kernel;
40
41namespace Mantid::Crystal {
42
43DECLARE_ALGORITHM(SCDCalibratePanels)
44
45// Default constructor
47 useAlgorithm("SCDCalibratePanels", 2);
48 deprecatedDate("2025-05-27");
49}
50
51const std::string SCDCalibratePanels::name() const { return "SCDCalibratePanels"; }
52
53int SCDCalibratePanels::version() const { return 1; }
54
55const std::string SCDCalibratePanels::category() const { return "Crystal\\Corrections"; }
56
58 PeaksWorkspace_sptr peaksWs = getProperty("PeakWorkspace");
59 // We must sort the peaks
60 std::vector<std::pair<std::string, bool>> criteria{{"BankName", true}};
61 peaksWs->sort(criteria);
62 // Remove peaks on edge
63 int edge = this->getProperty("EdgePixels");
64 Geometry::Instrument_const_sptr inst = peaksWs->getInstrument();
65 Geometry::ComponentInfo const &compInfo = peaksWs->componentInfo();
66 if (edge > 0) {
67 std::vector<Peak> &peaks = peaksWs->getPeaks();
68 auto it = std::remove_if(peaks.begin(), peaks.end(), [edge, &compInfo](const Peak &pk) {
69 return edgePixel(compInfo, pk.getBankName(), pk.getCol(), pk.getRow(), edge);
70 });
71 peaks.erase(it, peaks.end());
72 }
73 findU(peaksWs);
74
75 auto nPeaks = static_cast<int>(peaksWs->getNumberPeaks());
76 bool changeL1 = getProperty("ChangeL1");
77 bool changeT0 = getProperty("ChangeT0");
78 bool bankPanels = getProperty("CalibrateBanks");
79 bool snapPanels = getProperty("CalibrateSNAPPanels");
80
81 if (changeT0)
82 findT0(nPeaks, peaksWs);
83 if (changeL1)
84 findL1(nPeaks, peaksWs);
85
86 boost::container::flat_set<string> MyBankNames;
87 boost::container::flat_set<string> MyPanels;
88 if (snapPanels) {
89 MyPanels.insert("East");
90 MyPanels.insert("West");
91 int maxRecurseDepth = 4;
92
93 PRAGMA_OMP(parallel for schedule(dynamic, 1) )
94 for (int num = 1; num < 64; ++num) {
96 std::ostringstream mess;
97 mess << "bank" << num;
98 IComponent_const_sptr comp = inst->getComponentByName(mess.str(), maxRecurseDepth);
99 PARALLEL_CRITICAL(MyBankNames)
100 if (comp)
101 MyBankNames.insert(mess.str());
103 }
105 } else {
106 for (int i = 0; i < nPeaks; ++i) {
107 std::string bankName = peaksWs->getPeak(i).getBankName();
108 if (bankName != "None")
109 MyBankNames.insert(bankName);
110 }
111 }
112
113 std::vector<std::string> fit_workspaces(MyBankNames.size() + MyPanels.size(), "fit_");
114 std::vector<std::string> parameter_workspaces(MyBankNames.size() + MyPanels.size(), "params_");
115 int bankAndPanelCount = 0;
116 for (const auto &MyPanel : MyPanels) {
117 fit_workspaces[bankAndPanelCount] += MyPanel;
118 parameter_workspaces[bankAndPanelCount] += MyPanel;
119 bankAndPanelCount++;
120 }
121 if (snapPanels) {
122 findL2(MyPanels, peaksWs);
123 ITableWorkspace_sptr results = AnalysisDataService::Instance().retrieveWS<ITableWorkspace>("params_West");
124 double delta = results->cell<double>(4, 1);
125 g_log.notice() << "For west rotation change det_arc1 " << delta << " degrees\n";
126 results = AnalysisDataService::Instance().retrieveWS<ITableWorkspace>("params_East");
127 delta = results->cell<double>(4, 1);
128 g_log.notice() << "For east rotation change det_arc2 " << delta << " degrees\n";
129 }
130
131 for (const auto &MyBankName : MyBankNames) {
132 fit_workspaces[bankAndPanelCount] += MyBankName;
133 parameter_workspaces[bankAndPanelCount] += MyBankName;
134 bankAndPanelCount++;
135 }
136 if (bankPanels) {
137 findL2(MyBankNames, peaksWs);
138 }
139
140 // remove skipped banks
141 for (int j = bankAndPanelCount - 1; j >= 0; j--) {
142 if (!AnalysisDataService::Instance().doesExist(fit_workspaces[j]))
143 fit_workspaces.erase(fit_workspaces.begin() + j);
144 if (!AnalysisDataService::Instance().doesExist(parameter_workspaces[j]))
145 parameter_workspaces.erase(parameter_workspaces.begin() + j);
146 }
147
148 // Try again to optimize L1
149 if (changeL1) {
150 findL1(nPeaks, peaksWs);
151 parameter_workspaces.emplace_back("params_L1");
152 fit_workspaces.emplace_back("fit_L1");
153 }
154 // Add T0 files to groups
155 if (changeT0) {
156 parameter_workspaces.emplace_back("params_T0");
157 fit_workspaces.emplace_back("fit_T0");
158 }
159 std::sort(parameter_workspaces.begin(), parameter_workspaces.end());
160 std::sort(fit_workspaces.begin(), fit_workspaces.end());
161
162 // collect output of fit for each spectrum into workspace groups
163 auto groupAlg = AlgorithmManager::Instance().createUnmanaged("GroupWorkspaces");
164 groupAlg->initialize();
165 groupAlg->setProperty("InputWorkspaces", parameter_workspaces);
166 groupAlg->setProperty("OutputWorkspace", "Fit_Parameters");
167 groupAlg->execute();
168
169 groupAlg = AlgorithmManager::Instance().createUnmanaged("GroupWorkspaces");
170 groupAlg->initialize();
171 groupAlg->setProperty("InputWorkspaces", fit_workspaces);
172 groupAlg->setProperty("OutputWorkspace", "Fit_Residuals");
173 groupAlg->execute();
174
175 // Use new instrument for PeaksWorkspace
176 Geometry::Instrument_sptr inst2 = std::const_pointer_cast<Geometry::Instrument>(peaksWs->getInstrument());
177 Geometry::OrientedLattice lattice0 = peaksWs->mutableSample().getOrientedLattice();
179 for (int i = 0; i < nPeaks; i++) {
181 DataObjects::Peak &peak = peaksWs->getPeak(i);
182 try {
183 peak.setInstrument(inst2);
184 } catch (const std::exception &exc) {
185 g_log.notice() << "Problem in applying calibration to peak " << i << " : " << exc.what() << "\n";
186 }
188 }
190
191 // Find U again for optimized geometry and index peaks
192 findU(peaksWs);
193 // Save as DetCal and XML if requested
194 string DetCalFileName = getProperty("DetCalFilename");
195 API::Run &run = peaksWs->mutableRun();
196 double T0 = 0.0;
197 if (run.hasProperty("T0")) {
198 T0 = run.getPropertyValueAsType<double>("T0");
199 }
200 saveIsawDetCal(inst2, MyBankNames, T0, DetCalFileName);
201 string XmlFileName = getProperty("XmlFilename");
202 saveXmlFile(XmlFileName, MyBankNames, *inst2);
203 // create table of theoretical vs calculated
204 //----------------- Calculate & Create Calculated vs Theoretical
205 // workspaces------------------,);
206 MatrixWorkspace_sptr ColWksp =
207 Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", MyBankNames.size(), nPeaks, nPeaks);
208 ColWksp->setInstrument(inst2);
209 MatrixWorkspace_sptr RowWksp =
210 Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", MyBankNames.size(), nPeaks, nPeaks);
211 RowWksp->setInstrument(inst2);
212 MatrixWorkspace_sptr TofWksp =
213 Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", MyBankNames.size(), nPeaks, nPeaks);
214 TofWksp->setInstrument(inst2);
215 OrientedLattice lattice = peaksWs->mutableSample().getOrientedLattice();
216 const DblMatrix &UB = lattice.getUB();
217 // sort again since edge peaks can trace to other banks
218 peaksWs->sort(criteria);
219 PARALLEL_FOR_IF(Kernel::threadSafe(*ColWksp, *RowWksp, *TofWksp))
220 for (int i = 0; i < static_cast<int>(MyBankNames.size()); ++i) {
222 const std::string &bankName = *std::next(MyBankNames.begin(), i);
223 size_t k = bankName.find_last_not_of("0123456789");
224 int bank = 0;
225 if (k < bankName.length())
226 bank = boost::lexical_cast<int>(bankName.substr(k + 1));
227 ColWksp->getSpectrum(i).setSpectrumNo(specnum_t(bank));
228 RowWksp->getSpectrum(i).setSpectrumNo(specnum_t(bank));
229 TofWksp->getSpectrum(i).setSpectrumNo(specnum_t(bank));
230 auto &ColX = ColWksp->mutableX(i);
231 auto &ColY = ColWksp->mutableY(i);
232 auto &RowX = RowWksp->mutableX(i);
233 auto &RowY = RowWksp->mutableY(i);
234 auto &TofX = TofWksp->mutableX(i);
235 auto &TofY = TofWksp->mutableY(i);
236 int icount = 0;
237 for (int j = 0; j < nPeaks; j++) {
238 Peak peak = peaksWs->getPeak(j);
239 if (peak.getBankName() == bankName) {
240 try {
241 V3D q_lab = (peak.getGoniometerMatrix() * UB) * peak.getHKL() * M_2_PI;
242 Peak theoretical(peak.getInstrument(), q_lab);
243 ColX[icount] = peak.getCol();
244 ColY[icount] = theoretical.getCol();
245 RowX[icount] = peak.getRow();
246 RowY[icount] = theoretical.getRow();
247 TofX[icount] = peak.getTOF();
248 TofY[icount] = theoretical.getTOF();
249 } catch (...) {
250 // g_log.debug() << "Problem only in printing peaks\n";
251 }
252 icount++;
253 }
254 }
256 }
258
259 string colFilename = getProperty("ColFilename");
260 string rowFilename = getProperty("RowFilename");
261 string tofFilename = getProperty("TofFilename");
262 saveNexus(colFilename, ColWksp);
263 saveNexus(rowFilename, RowWksp);
264 saveNexus(tofFilename, TofWksp);
265}
266
267void SCDCalibratePanels::saveNexus(const std::string &outputFile, const MatrixWorkspace_sptr &outputWS) {
268 auto save = createChildAlgorithm("SaveNexus");
269 save->setProperty("InputWorkspace", outputWS);
270 save->setProperty("FileName", outputFile);
271 save->execute();
272}
273
275 MatrixWorkspace_sptr L1WS = std::dynamic_pointer_cast<MatrixWorkspace>(
276 API::WorkspaceFactory::Instance().create("Workspace2D", 1, 3 * nPeaks, 3 * nPeaks));
277
278 IAlgorithm_sptr fitL1_alg;
279 try {
280 fitL1_alg = createChildAlgorithm("Fit", -1, -1, false);
281 } catch (Exception::NotFoundError &) {
282 g_log.error("Can't locate Fit algorithm");
283 throw;
284 }
285 std::ostringstream fun_str;
286 fun_str << "name=SCDPanelErrors,Workspace=" << peaksWs->getName() << ",Bank=moderator";
287 std::ostringstream tie_str;
288 tie_str << "XShift=0.0,YShift=0.0,XRotate=0.0,YRotate=0.0,ZRotate=0.0,"
289 "ScaleWidth=1.0,ScaleHeight=1.0,T0Shift ="
290 << mT0;
291 fitL1_alg->setPropertyValue("Function", fun_str.str());
292 fitL1_alg->setProperty("Ties", tie_str.str());
293 fitL1_alg->setProperty("InputWorkspace", L1WS);
294 fitL1_alg->setProperty("CreateOutput", true);
295 fitL1_alg->setProperty("Output", "fit");
296 fitL1_alg->executeAsChildAlg();
297 std::string fitL1Status = fitL1_alg->getProperty("OutputStatus");
298 double chisqL1 = fitL1_alg->getProperty("OutputChi2overDoF");
299 MatrixWorkspace_sptr fitL1 = fitL1_alg->getProperty("OutputWorkspace");
300 AnalysisDataService::Instance().addOrReplace("fit_L1", fitL1);
301 ITableWorkspace_sptr paramsL1 = fitL1_alg->getProperty("OutputParameters");
302 AnalysisDataService::Instance().addOrReplace("params_L1", paramsL1);
303 double deltaL1 = paramsL1->getRef<double>("Value", 2);
304 SCDPanelErrors com;
305 com.moveDetector(0.0, 0.0, deltaL1, 0.0, 0.0, 0.0, 1.0, 1.0, "moderator", peaksWs);
306 g_log.notice() << "L1 = " << -peaksWs->getInstrument()->getSource()->getPos().Z() << " " << fitL1Status
307 << " Chi2overDoF " << chisqL1 << "\n";
308}
309
311 MatrixWorkspace_sptr T0WS = std::dynamic_pointer_cast<MatrixWorkspace>(
312 API::WorkspaceFactory::Instance().create("Workspace2D", 1, 3 * nPeaks, 3 * nPeaks));
313
314 IAlgorithm_sptr fitT0_alg;
315 try {
316 fitT0_alg = createChildAlgorithm("Fit", -1, -1, false);
317 } catch (Exception::NotFoundError &) {
318 g_log.error("Can't locate Fit algorithm");
319 throw;
320 }
321 std::ostringstream fun_str;
322 fun_str << "name=SCDPanelErrors,Workspace=" << peaksWs->getName() << ",Bank=none";
323 std::ostringstream tie_str;
324 tie_str << "XShift=0.0,YShift=0.0,ZShift=0.0,XRotate=0.0,YRotate=0.0,ZRotate=0.0,"
325 "ScaleWidth=1.0,ScaleHeight=1.0";
326 fitT0_alg->setPropertyValue("Function", fun_str.str());
327 fitT0_alg->setProperty("Ties", tie_str.str());
328 fitT0_alg->setProperty("InputWorkspace", T0WS);
329 fitT0_alg->setProperty("CreateOutput", true);
330 fitT0_alg->setProperty("Output", "fit");
331 // Does not converge with derviative minimizers
332 fitT0_alg->setProperty("Minimizer", "Simplex");
333 fitT0_alg->setProperty("MaxIterations", 1000);
334 fitT0_alg->executeAsChildAlg();
335 std::string fitT0Status = fitT0_alg->getProperty("OutputStatus");
336 double chisqT0 = fitT0_alg->getProperty("OutputChi2overDoF");
337 MatrixWorkspace_sptr fitT0 = fitT0_alg->getProperty("OutputWorkspace");
338 AnalysisDataService::Instance().addOrReplace("fit_T0", fitT0);
339 ITableWorkspace_sptr paramsT0 = fitT0_alg->getProperty("OutputParameters");
340 AnalysisDataService::Instance().addOrReplace("params_T0", paramsT0);
341 mT0 = paramsT0->getRef<double>("Value", 8);
342 API::Run &run = peaksWs->mutableRun();
343 // set T0 in the run parameters adding to value in peaks file
344 double oldT0 = 0.0;
345 if (run.hasProperty("T0")) {
346 oldT0 = run.getPropertyValueAsType<double>("T0");
347 }
348 run.addProperty<double>("T0", mT0 + oldT0, true);
349 g_log.notice() << "T0 = " << mT0 << " " << fitT0Status << " Chi2overDoF " << chisqT0 << "\n";
350 for (int i = 0; i < peaksWs->getNumberPeaks(); i++) {
351 DataObjects::Peak &peak = peaksWs->getPeak(i);
352
354
355 wl.initialize(peak.getL1(), 0, {{UnitParams::l2, peak.getL2()}, {UnitParams::twoTheta, peak.getScattering()}});
356 peak.setWavelength(wl.singleFromTOF(peak.getTOF() + mT0));
357 }
358}
359
361 IAlgorithm_sptr ub_alg;
362 try {
363 ub_alg = createChildAlgorithm("CalculateUMatrix", -1, -1, false);
364 } catch (Exception::NotFoundError &) {
365 g_log.error("Can't locate CalculateUMatrix algorithm");
366 throw;
367 }
368 double a = getProperty("a");
369 double b = getProperty("b");
370 double c = getProperty("c");
371 double alpha = getProperty("alpha");
372 double beta = getProperty("beta");
373 double gamma = getProperty("gamma");
374 if ((a == EMPTY_DBL() || b == EMPTY_DBL() || c == EMPTY_DBL() || alpha == EMPTY_DBL() || beta == EMPTY_DBL() ||
375 gamma == EMPTY_DBL()) &&
376 peaksWs->sample().hasOrientedLattice()) {
377 OrientedLattice latt = peaksWs->mutableSample().getOrientedLattice();
378 a = latt.a();
379 b = latt.b();
380 c = latt.c();
381 alpha = latt.alpha();
382 beta = latt.beta();
383 gamma = latt.gamma();
384 }
385 ub_alg->setProperty("PeaksWorkspace", peaksWs);
386 ub_alg->setProperty("a", a);
387 ub_alg->setProperty("b", b);
388 ub_alg->setProperty("c", c);
389 ub_alg->setProperty("alpha", alpha);
390 ub_alg->setProperty("beta", beta);
391 ub_alg->setProperty("gamma", gamma);
392 ub_alg->executeAsChildAlg();
393
394 // Reindex peaks with new UB
395 auto alg = createChildAlgorithm("IndexPeaks");
396 alg->setPropertyValue("PeaksWorkspace", peaksWs->getName());
397 alg->setProperty("Tolerance", 0.15);
398 alg->executeAsChildAlg();
399 int numIndexed = alg->getProperty("NumIndexed");
400 g_log.notice() << "Number Indexed = " << numIndexed << "\n";
401 g_log.notice() << peaksWs->sample().getOrientedLattice().getUB() << "\n";
402}
403
414void SCDCalibratePanels::saveIsawDetCal(std::shared_ptr<Instrument> &instrument,
415 boost::container::flat_set<string> &AllBankName, double T0,
416 const string &filename) {
417 // having a filename triggers doing the work
418 if (filename.empty())
419 return;
420
421 g_log.notice() << "Saving DetCal file in " << filename << "\n";
422
423 // create a workspace to pass to SaveIsawDetCal
424 const size_t number_spectra = instrument->getNumberDetectors();
425 DataObjects::Workspace2D_sptr wksp = std::dynamic_pointer_cast<DataObjects::Workspace2D>(
426 WorkspaceFactory::Instance().create("Workspace2D", number_spectra, 2, 1));
427 wksp->setInstrument(instrument);
428 wksp->rebuildSpectraMapping(true /* include monitors */);
429
430 // convert the bank names into a vector
431 std::vector<string> banknames(AllBankName.begin(), AllBankName.end());
432
433 // call SaveIsawDetCal
434 auto alg = createChildAlgorithm("SaveIsawDetCal");
435 alg->setProperty("InputWorkspace", wksp);
436 alg->setProperty("Filename", filename);
437 alg->setProperty("TimeOffset", T0);
438 alg->setProperty("BankNames", banknames);
439 alg->executeAsChildAlg();
440}
441
443 declareProperty(std::make_unique<WorkspaceProperty<PeaksWorkspace>>("PeakWorkspace", "", Kernel::Direction::InOut),
444 "Workspace of Indexed Peaks");
445
446 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
447 mustBePositive->setLower(0.0);
448
449 declareProperty("a", EMPTY_DBL(), mustBePositive,
450 "Lattice Parameter a (Leave empty to use lattice constants "
451 "in peaks workspace)");
452 declareProperty("b", EMPTY_DBL(), mustBePositive,
453 "Lattice Parameter b (Leave empty to use lattice constants "
454 "in peaks workspace)");
455 declareProperty("c", EMPTY_DBL(), mustBePositive,
456 "Lattice Parameter c (Leave empty to use lattice constants "
457 "in peaks workspace)");
458 declareProperty("alpha", EMPTY_DBL(), mustBePositive,
459 "Lattice Parameter alpha in degrees (Leave empty to use "
460 "lattice constants in peaks workspace)");
461 declareProperty("beta", EMPTY_DBL(), mustBePositive,
462 "Lattice Parameter beta in degrees (Leave empty to use "
463 "lattice constants in peaks workspace)");
464 declareProperty("gamma", EMPTY_DBL(), mustBePositive,
465 "Lattice Parameter gamma in degrees (Leave empty to use "
466 "lattice constants in peaks workspace)");
467 declareProperty("ChangeL1", true, "Change the L1(source to sample) distance");
468 declareProperty("ChangeT0", false, "Change the T0 (initial TOF)");
469 declareProperty("ChangePanelSize", true,
470 "Change the height and width of the "
471 "detectors. Implemented only for "
472 "RectangularDetectors.");
473
474 declareProperty("EdgePixels", 0, "Remove peaks that are at pixels this close to edge. ");
475 declareProperty("CalibrateBanks", true, "Calibrate the panels of the banks.");
476 declareProperty("CalibrateSNAPPanels", false,
477 "Calibrate the 3 X 3 panels of the "
478 "sides of SNAP.");
479
480 // ---------- outputs
481 const std::vector<std::string> detcalExts{".DetCal", ".Det_Cal"};
483 std::make_unique<FileProperty>("DetCalFilename", "SCDCalibrate.DetCal", FileProperty::Save, detcalExts),
484 "Path to an ISAW-style .detcal file to save.");
485
486 declareProperty(std::make_unique<FileProperty>("XmlFilename", "", FileProperty::OptionalSave, ".xml"),
487 "Path to an Mantid .xml description(for LoadParameterFile) file to "
488 "save.");
489
490 declareProperty(std::make_unique<FileProperty>("ColFilename", "ColCalcvsTheor.nxs", FileProperty::Save, ".nxs"),
491 "Path to a NeXus file comparing calculated and theoretical "
492 "column of each peak.");
493
494 declareProperty(std::make_unique<FileProperty>("RowFilename", "RowCalcvsTheor.nxs", FileProperty::Save, ".nxs"),
495 "Path to a NeXus file comparing calculated and theoretical "
496 "row of each peak.");
497
498 declareProperty(std::make_unique<FileProperty>("TofFilename", "TofCalcvsTheor.nxs", FileProperty::Save, ".nxs"),
499 "Path to a NeXus file comparing calculated and theoretical "
500 "TOF of each peak.");
501
502 const string OUTPUTS("Outputs");
503 setPropertyGroup("DetCalFilename", OUTPUTS);
504 setPropertyGroup("XmlFilename", OUTPUTS);
505 setPropertyGroup("ColFilename", OUTPUTS);
506 setPropertyGroup("RowFilename", OUTPUTS);
507 setPropertyGroup("TofFilename", OUTPUTS);
508}
509
510void writeXmlParameter(ofstream &ostream, const string &name, const double value) {
511 ostream << " <parameter name =\"" << name << "\"><value val=\"" << value << "\" /> </parameter>\n";
512}
513
514void SCDCalibratePanels::saveXmlFile(const string &FileName, const boost::container::flat_set<string> &AllBankNames,
515 const Instrument &instrument) const {
516 if (FileName.empty())
517 return;
518
519 g_log.notice() << "Saving parameter file as " << FileName << "\n";
520
521 // create the file and add the header
522 ofstream oss3(FileName.c_str());
523 oss3 << "<?xml version=\"1.0\" encoding=\"UTF-8\" ?>\n";
524 oss3 << " <parameter-file instrument=\"" << instrument.getName() << "\" valid-from=\""
525 << instrument.getValidFromDate().toISO8601String() << "\">\n";
527
528 // write out the detector banks
529 for (auto bankName : AllBankNames) {
530 if (instrument.getName().compare("CORELLI") == 0.0)
531 bankName.append("/sixteenpack");
532 oss3 << "<component-link name=\"" << bankName << "\">\n";
533 std::shared_ptr<const IComponent> bank = instrument.getComponentByName(bankName);
534
535 Quat relRot = bank->getRelativeRot();
536
537 std::vector<double> relRotAngles = relRot.getEulerAngles("XYZ");
538
539 writeXmlParameter(oss3, "rotx", relRotAngles[0]);
540 writeXmlParameter(oss3, "roty", relRotAngles[1]);
541 writeXmlParameter(oss3, "rotz", relRotAngles[2]);
542
543 V3D pos1 = bank->getRelativePos();
544 writeXmlParameter(oss3, "x", pos1.X());
545 writeXmlParameter(oss3, "y", pos1.Y());
546 writeXmlParameter(oss3, "z", pos1.Z());
547
548 vector<double> oldScalex = pmap->getDouble(bank->getName(), string("scalex"));
549 vector<double> oldScaley = pmap->getDouble(bank->getName(), string("scaley"));
550
551 double scalex, scaley;
552 if (!oldScalex.empty())
553 scalex = oldScalex[0];
554 else
555 scalex = 1.;
556
557 if (!oldScaley.empty())
558 scaley = oldScaley[0];
559 else
560 scaley = 1.;
561
562 oss3 << R"( <parameter name ="scalex"><value val=")" << scalex << "\" /> </parameter>\n";
563 oss3 << R"( <parameter name ="scaley"><value val=")" << scaley << "\" /> </parameter>\n";
564 oss3 << "</component-link>\n";
565 } // for each bank in the group
566
567 // write out the source
568 IComponent_const_sptr source = instrument.getSource();
569
570 oss3 << "<component-link name=\"" << source->getName() << "\">\n";
571 V3D sourceRelPos = source->getRelativePos();
572
573 writeXmlParameter(oss3, "x", sourceRelPos.X());
574 writeXmlParameter(oss3, "y", sourceRelPos.Y());
575 writeXmlParameter(oss3, "z", sourceRelPos.Z());
576 oss3 << "</component-link>\n";
577 oss3 << "</parameter-file>\n";
578
579 // flush and close the file
580 oss3.flush();
581 oss3.close();
582}
583void SCDCalibratePanels::findL2(boost::container::flat_set<string> MyBankNames,
584 const DataObjects::PeaksWorkspace_sptr &peaksWs) {
585 bool changeSize = getProperty("ChangePanelSize");
586 Geometry::Instrument_const_sptr inst = peaksWs->getInstrument();
587
589 for (int bankIndex = 0; bankIndex < static_cast<int>(MyBankNames.size()); ++bankIndex) {
591 const std::string &iBank = *std::next(MyBankNames.begin(), bankIndex);
592 const std::string bankName = "__PWS_" + iBank;
593 PeaksWorkspace_sptr local = peaksWs->clone();
594 AnalysisDataService::Instance().addOrReplace(bankName, local);
595 std::vector<Peak> &localPeaks = local->getPeaks();
596 auto lit = std::remove_if(localPeaks.begin(), localPeaks.end(), [&iBank](const Peak &pk) {
597 std::string name = pk.getBankName();
598 IComponent_const_sptr det = pk.getInstrument()->getComponentByName(name);
599 if (det && iBank.substr(0, 4) != "bank") {
600 IComponent_const_sptr parent = det->getParent();
601 if (parent) {
602 IComponent_const_sptr grandparent = parent->getParent();
603 if (grandparent) {
604 name = grandparent->getName();
605 }
606 }
607 }
608
609 return name != iBank;
610 });
611 localPeaks.erase(lit, localPeaks.end());
612
613 int nBankPeaks = local->getNumberPeaks();
614 if (nBankPeaks < 6) {
615 g_log.notice() << "Too few peaks for " << iBank << "\n";
616 continue;
617 }
618
619 MatrixWorkspace_sptr q3DWS = std::dynamic_pointer_cast<MatrixWorkspace>(
620 API::WorkspaceFactory::Instance().create("Workspace2D", 1, 3 * nBankPeaks, 3 * nBankPeaks));
621
622 auto &outSpec = q3DWS->getSpectrum(0);
623 auto &yVec = outSpec.mutableY();
624 auto &eVec = outSpec.mutableE();
625 auto &xVec = outSpec.mutableX();
626 yVec = 0.0;
627
628 for (int i = 0; i < nBankPeaks; i++) {
629 const DataObjects::Peak &peak = local->getPeak(i);
630 // 1/sigma is considered the weight for the fit
631 double weight = 1.; // default is even weighting
632 if (peak.getSigmaIntensity() > 0.) // prefer weight by sigmaI
633 weight = 1.0 / peak.getSigmaIntensity();
634 else if (peak.getIntensity() > 0.) // next favorite weight by I
635 weight = 1.0 / peak.getIntensity();
636 else if (peak.getBinCount() > 0.) // then by counts in peak centre
637 weight = 1.0 / peak.getBinCount();
638 for (int j = 0; j < 3; j++) {
639 int k = i * 3 + j;
640 xVec[k] = k;
641 eVec[k] = weight;
642 }
643 }
644
645 IAlgorithm_sptr fit_alg;
646 try {
647 fit_alg = createChildAlgorithm("Fit", -1, -1, false);
648 } catch (Exception::NotFoundError &) {
649 g_log.error("Can't locate Fit algorithm");
650 throw;
651 }
652 std::ostringstream fun_str;
653 fun_str << "name=SCDPanelErrors,Workspace=" + bankName << ",Bank=" << iBank;
654 fit_alg->setPropertyValue("Function", fun_str.str());
655 std::ostringstream tie_str;
656 tie_str << "ScaleWidth=1.0,ScaleHeight=1.0,T0Shift =" << mT0;
657 fit_alg->setProperty("Ties", tie_str.str());
658 fit_alg->setProperty("InputWorkspace", q3DWS);
659 fit_alg->setProperty("CreateOutput", true);
660 fit_alg->setProperty("Output", "fit");
661 fit_alg->executeAsChildAlg();
662 std::string fitStatus = fit_alg->getProperty("OutputStatus");
663 double fitChisq = fit_alg->getProperty("OutputChi2overDoF");
664 g_log.notice() << iBank << " " << fitStatus << " Chi2overDoF " << fitChisq << "\n";
665 MatrixWorkspace_sptr fitWS = fit_alg->getProperty("OutputWorkspace");
666 AnalysisDataService::Instance().addOrReplace("fit_" + iBank, fitWS);
667 ITableWorkspace_sptr paramsWS = fit_alg->getProperty("OutputParameters");
668 AnalysisDataService::Instance().addOrReplace("params_" + iBank, paramsWS);
669 double xShift = paramsWS->getRef<double>("Value", 0);
670 double yShift = paramsWS->getRef<double>("Value", 1);
671 double zShift = paramsWS->getRef<double>("Value", 2);
672 double xRotate = paramsWS->getRef<double>("Value", 3);
673 double yRotate = paramsWS->getRef<double>("Value", 4);
674 double zRotate = paramsWS->getRef<double>("Value", 5);
675 double scaleWidth = 1.0;
676 double scaleHeight = 1.0;
677 // Scaling only implemented for Rectangular Detectors
678 Geometry::IComponent_const_sptr comp = peaksWs->getInstrument()->getComponentByName(iBank);
679 std::shared_ptr<const Geometry::RectangularDetector> rectDet =
680 std::dynamic_pointer_cast<const Geometry::RectangularDetector>(comp);
681 if (rectDet && changeSize) {
682 IAlgorithm_sptr fit2_alg;
683 try {
684 fit2_alg = createChildAlgorithm("Fit", -1, -1, false);
685 } catch (Exception::NotFoundError &) {
686 g_log.error("Can't locate Fit algorithm");
687 throw;
688 }
689 fit2_alg->setPropertyValue("Function", fun_str.str());
690 std::ostringstream tie_str2;
691 tie_str2 << "XShift=" << xShift << ",YShift=" << yShift << ",ZShift=" << zShift << ",XRotate=" << xRotate
692 << ",YRotate=" << yRotate << ",ZRotate=" << zRotate << ",T0Shift =" << mT0;
693 fit2_alg->setProperty("Ties", tie_str2.str());
694 fit2_alg->setProperty("InputWorkspace", q3DWS);
695 fit2_alg->setProperty("CreateOutput", true);
696 fit2_alg->setProperty("Output", "fit");
697 fit2_alg->executeAsChildAlg();
698 std::string fit2Status = fit2_alg->getProperty("OutputStatus");
699 double fit2Chisq = fit2_alg->getProperty("OutputChi2overDoF");
700 g_log.notice() << iBank << " " << fit2Status << " Chi2overDoF " << fit2Chisq << "\n";
701 fitWS = fit2_alg->getProperty("OutputWorkspace");
702 AnalysisDataService::Instance().addOrReplace("fit_" + iBank, fitWS);
703 paramsWS = fit2_alg->getProperty("OutputParameters");
704 AnalysisDataService::Instance().addOrReplace("params_" + iBank, paramsWS);
705 scaleWidth = paramsWS->getRef<double>("Value", 6);
706 scaleHeight = paramsWS->getRef<double>("Value", 7);
707 }
708 AnalysisDataService::Instance().remove(bankName);
709 SCDPanelErrors det;
710 det.moveDetector(xShift, yShift, zShift, xRotate, yRotate, zRotate, scaleWidth, scaleHeight, iBank, peaksWs);
712 }
714}
715} // namespace Mantid::Crystal
std::string name
Definition Run.cpp:60
#define DECLARE_ALGORITHM(classname)
Definition Algorithm.h:538
double value
The value of the point.
Definition FitMW.cpp:51
int numIndexed
#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...
#define PARALLEL_CRITICAL(name)
#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 PRAGMA_OMP(expression)
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
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.
Kernel::Logger & g_log
Definition Algorithm.h:422
@ OptionalSave
to specify a file to write to but an empty string is
@ Save
to specify a file to write to, the file may or may not exist
ITableWorkspace is an implementation of Workspace in which the data are organised in columns of same ...
T & cell(size_t row, size_t col)
Get the reference to the element in row row and column col.
bool hasProperty(const std::string &name) const
Does the property exist on the object.
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
This class stores information regarding an experimental run as a series of log entries.
Definition Run.h:35
A property class for workspaces.
SCDCalibratePanels calibrates instrument parameters for Rectangular Detectors.
void saveXmlFile(const std::string &FileName, const boost::container::flat_set< std::string > &AllBankNames, const Geometry::Instrument &instrument) const
Saves the new instrument to an xml file that can be used with the LoadParameterFile Algorithm.
void init() override
Virtual method - must be overridden by concrete algorithm.
const std::string name() const override
function to return a name of the algorithm, must be overridden in all algorithms
void findT0(int nPeaks, const DataObjects::PeaksWorkspace_sptr &peaksWs)
Function to optimize T0.
void findL2(boost::container::flat_set< std::string > MyBankNames, const DataObjects::PeaksWorkspace_sptr &peaksWs)
Function to optimize L2.
int version() const override
Algorithm's version for identification overriding a virtual method.
void saveIsawDetCal(std::shared_ptr< Geometry::Instrument > &instrument, boost::container::flat_set< std::string > &AllBankName, double T0, const std::string &filename)
Really this is the operator SaveIsawDetCal but only the results of the given banks are saved.
void findU(const DataObjects::PeaksWorkspace_sptr &peaksWs)
Function to calculate U.
void exec() override
Virtual method - must be overridden by concrete algorithm.
void saveNexus(const std::string &outputFile, const API::MatrixWorkspace_sptr &outputWS)
save workspaces
void findL1(int nPeaks, const DataObjects::PeaksWorkspace_sptr &peaksWs)
Function to optimize L1.
const std::string category() const override
Algorithm's category for identification overriding a virtual method.
void moveDetector(double x, double y, double z, double rotx, double roty, double rotz, double scalex, double scaley, std::string detname, const API::Workspace_sptr &inputW) const
Move detectors with parameters.
double getIntensity() const override
Return the integrated peak intensity.
Definition BasePeak.cpp:184
double getSigmaIntensity() const override
Return the error on the integrated peak intensity.
Definition BasePeak.cpp:187
Mantid::Kernel::V3D getHKL() const override
Return the HKL vector.
Definition BasePeak.cpp:102
double getBinCount() const override
Return the # of counts in the bin at its peak.
Definition BasePeak.cpp:181
Mantid::Kernel::Matrix< double > getGoniometerMatrix() const override
Get the goniometer rotation matrix at which this peak was measured.
Definition BasePeak.cpp:209
Structure describing a single-crystal peak.
Definition Peak.h:34
int getCol() const override
For RectangularDetectors only, returns the column (x) of the pixel of the detector or -1 if not found...
Definition Peak.cpp:335
Geometry::Instrument_const_sptr getInstrument() const
Return a shared ptr to the instrument for this peak.
Definition Peak.cpp:320
int getRow() const override
For RectangularDetectors only, returns the row (y) of the pixel of the detector or -1 if not found.
Definition Peak.cpp:326
void setWavelength(double wavelength) override
Set the incident wavelength of the neutron.
Definition Peak.cpp:188
double getL1() const override
Return the L1 flight path length (source to sample), in meters.
Definition Peak.cpp:724
double getTOF() const override
Calculate the time of flight (in microseconds) of the neutrons for this peak, using the geometry of t...
Definition Peak.cpp:373
double getL2() const override
Return the L2 flight path length (sample to detector), in meters.
Definition Peak.cpp:728
const std::string & getBankName() const
Find the name of the bank that is the parent of the detector.
Definition Peak.cpp:347
void setInstrument(const Geometry::Instrument_const_sptr &inst)
Set the instrument (and save the source/sample pos).
Definition Peak.cpp:294
double getScattering() const override
Calculate the scattering angle of the peak
Definition Peak.cpp:391
std::shared_ptr< const IComponent > getComponentByName(const std::string &cname, int nlevels=0) const override
Returns a pointer to the first component of assembly encountered with the given name.
ComponentInfo : Provides a component centric view on to the instrument.
std::string getName() const override
Get the IComponent name.
Base Instrument Class.
Definition Instrument.h:48
IComponent_const_sptr getSource() const
Gets a pointer to the source.
std::shared_ptr< ParameterMap > getParameterMap() const
Pointer to the NOT const ParameterMap holding the parameters of the modified instrument components.
Types::Core::DateAndTime getValidFromDate() const
Definition Instrument.h:171
Class to implement UB matrix.
const Kernel::DblMatrix & getUB() const
Get the 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
Exception for when an item is not found in a collection.
Definition Exception.h:145
void setPropertyGroup(const std::string &name, const std::string &group)
Set the group for a given property.
void notice(const std::string &msg)
Logs at notice level.
Definition Logger.cpp:126
void error(const std::string &msg)
Logs at error level.
Definition Logger.cpp:108
Class for quaternions.
Definition Quat.h:39
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...
void initialize(const double &_l1, const int &_emode, const UnitParametersMap &params)
Initialize the unit to perform conversion using singleToTof() and singleFromTof()
Definition Unit.cpp:133
Wavelength in Angstrom.
Definition Unit.h:267
double singleFromTOF(const double tof) const override
Convert a single tof value to this unit.
Definition Unit.cpp:394
Class for 3D vectors.
Definition V3D.h:34
constexpr double X() const noexcept
Get x.
Definition V3D.h:238
constexpr double Y() const noexcept
Get y.
Definition V3D.h:239
constexpr double Z() const noexcept
Get z.
Definition V3D.h:240
std::shared_ptr< IAlgorithm > IAlgorithm_sptr
shared pointer to Mantid::API::IAlgorithm
std::shared_ptr< ITableWorkspace > ITableWorkspace_sptr
shared pointer to Mantid::API::ITableWorkspace
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
void writeXmlParameter(ofstream &ostream, const string &name, const double value)
std::shared_ptr< PeaksWorkspace > PeaksWorkspace_sptr
Typedef for a shared pointer to a peaks workspace.
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:167
std::shared_ptr< ParameterMap > ParameterMap_sptr
ParameterMap shared pointer typedef.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
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.
int32_t specnum_t
Typedef for a spectrum Number.
Definition IDTypes.h:14
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Definition EmptyValues.h:42
Generate a tableworkspace to store the calibration results.
adjust instrument component position and orientation
: detector size scale at y-direction
STL namespace.
@ InOut
Both an input & output workspace.
Definition Property.h:55