From 3d564fe85f6865daf4582d35146e9352dbd0bf3f Mon Sep 17 00:00:00 2001
From: Vickie Lynch <lynchve@ornl.gov>
Date: Mon, 25 Jul 2016 18:52:09 -0400
Subject: [PATCH] Refs #16199 optimize one bank
---
.../Crystal/inc/MantidCrystal/PeakHKLErrors.h | 22 +
.../inc/MantidCrystal/SCDCalibratePanels.h | 2 +-
.../inc/MantidCrystal/SCDPanelErrors.h | 227 ++---
Framework/Crystal/src/PeakHKLErrors.cpp | 41 +-
Framework/Crystal/src/SCDCalibratePanels.cpp | 720 ++------------
Framework/Crystal/src/SCDPanelErrors.cpp | 905 +++++-------------
.../src/MoveInstrumentComponent.cpp | 46 +-
7 files changed, 474 insertions(+), 1489 deletions(-)
diff --git a/Framework/Crystal/inc/MantidCrystal/PeakHKLErrors.h b/Framework/Crystal/inc/MantidCrystal/PeakHKLErrors.h
index c5fb91c9fec..c9d3a65b398 100644
--- a/Framework/Crystal/inc/MantidCrystal/PeakHKLErrors.h
+++ b/Framework/Crystal/inc/MantidCrystal/PeakHKLErrors.h
@@ -17,6 +17,7 @@
#include "MantidAPI/IFunction.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/Matrix.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
namespace Mantid {
namespace Crystal {
@@ -65,6 +66,27 @@ public:
void init() override;
+ /**
+ * Creates a new peak, matching the old peak except for a different
+ *instrument.
+ *
+ * The Time of flightis the same except offset by T0. L0 should be the L0
+ *for the new instrument.
+ * It is added as a parameter in case the instrument will have the initial
+ *flight path adjusted later.
+ * NOTE: the wavelength is changed.
+ *
+ * @param peak_old - The old peak
+ * @param instrNew -The new instrument
+ * @param T0 :
+ * @param L0 :
+ * @return The new peak with the new instrument( adjusted with the
+ *parameters) and time adjusted.
+ */
+ static DataObjects::Peak createNewPeak(const Geometry::IPeak &peak_old,
+ Geometry::Instrument_sptr instrNew,
+ double T0, double L0);
+
static void cLone(boost::shared_ptr<Geometry::ParameterMap> &pmap,
boost::shared_ptr<const Geometry::IComponent> component,
boost::shared_ptr<const Geometry::ParameterMap> &pmapSv);
diff --git a/Framework/Crystal/inc/MantidCrystal/SCDCalibratePanels.h b/Framework/Crystal/inc/MantidCrystal/SCDCalibratePanels.h
index d90e4605040..2f257bff0fc 100644
--- a/Framework/Crystal/inc/MantidCrystal/SCDCalibratePanels.h
+++ b/Framework/Crystal/inc/MantidCrystal/SCDCalibratePanels.h
@@ -238,7 +238,7 @@ public:
std::string bankPrefixName);
private:
- void saveIsawDetCal(boost::shared_ptr<const Geometry::Instrument> &instrument,
+ void saveIsawDetCal(boost::shared_ptr<Geometry::Instrument> &instrument,
std::set<std::string> &AllBankName, double T0,
std::string filename);
diff --git a/Framework/Crystal/inc/MantidCrystal/SCDPanelErrors.h b/Framework/Crystal/inc/MantidCrystal/SCDPanelErrors.h
index 60097b2e682..d4ac0d95406 100644
--- a/Framework/Crystal/inc/MantidCrystal/SCDPanelErrors.h
+++ b/Framework/Crystal/inc/MantidCrystal/SCDPanelErrors.h
@@ -4,201 +4,98 @@
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
-//#include "MantidCurveFitting/DllConfig.h"
#include "MantidAPI/ParamFunction.h"
#include "MantidAPI/IFunction1D.h"
-#include "MantidDataObjects/PeaksWorkspace.h"
-#include "MantidDataObjects/Workspace2D.h"
-#include "MantidGeometry/Crystal/UnitCell.h"
-#include <boost/lexical_cast.hpp>
+#include "MantidKernel/System.h"
+#include "MantidAPI/Workspace_fwd.h"
+#include <cmath>
namespace Mantid {
namespace Crystal {
+/*
+@author Vickie Lynch, SNS
+@date 7/25/2016
-/**
- * Copyright © 2011-12 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
- *National Laboratory & European Spallation Source
- *
- * This file is part of Mantid.
- *
- * Mantid is free software; you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation; either version 3 of the License, or
- * (at your option) any later version.
- *
- * Mantid is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- * File change history is stored at: <https://github.com/mantidproject/mantid>
- * Code Documentation is available at: <http://doxygen.mantidproject.org>
- */
+Copyright © 2007-8 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+National Laboratory & European Spallation Source
+This file is part of Mantid.
+
+Mantid is free software; you can redistribute it and/or modify
+it under the terms of the GNU General Public License as published by
+the Free Software Foundation; either version 3 of the License, or
+(at your option) any later version.
+
+Mantid is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License
+along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+File change history is stored at: <https://github.com/mantidproject/mantid>
+Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
class DLLExport SCDPanelErrors : public API::ParamFunction,
- public API::IFunction1D {
+ public API::IFunction1D {
public:
+ /// Constructor
SCDPanelErrors();
- /**
- * Constructor
- * @param pwk - The PeaksWorkspace
- * @param Component_name - The comma separated list of bank names for which
- *this
- *Function
- * this Function calculates the associated errors in
- *qx,qy,qz
- * @param ax - The a lattice parameter
- * @param bx - The b lattice parameter
- * @param cx - The c lattice parameter
- * @param alphax - The alpha lattice parameter in degrees
- * @param betax - The beta lattice parameter in degrees
- * @param gammax - The gamma lattice parameter in degrees
- * @param tolerance1 - The maximum distance a peak's h value, k value and
- *lvalue are
- * from an integer to be considered indexed. Outside of
- *this
- * constructor, ALL PEAKS are considered INDEXED.
- *
- *
- */
- SCDPanelErrors(DataObjects::PeaksWorkspace_sptr &pwk,
- std::string &Component_name, double ax, double bx, double cx,
- double alphax, double betax, double gammax, double tolerance1);
-
+ /// overwrite IFunction base class methods
std::string name() const override { return "SCDPanelErrors"; }
-
- const std::string category() const override { return "Calibrate"; }
-
+ const std::string category() const override { return "General"; }
void function1D(double *out, const double *xValues,
const size_t nData) const override;
-
+ /// function derivatives
void functionDeriv1D(API::Jacobian *out, const double *xValues,
const size_t nData) override;
- size_t nAttributes() const override;
-
- std::vector<std::string> getAttributeNames() const override;
-
- Attribute getAttribute(const std::string &attName) const override;
-
+ /// Set a value to attribute attName
void setAttribute(const std::string &attName,
- const Attribute &value) override;
-
- bool hasAttribute(const std::string &attName) const override;
-
- /**
- * A utility method that will set up the Workspace needed by this function.
- * @param pwks -The peaksWorkspace. All peaks indexed to the given
- *tolerance and whose
- * associated bankName matches on of the strings in bankNames
- *will be included.
- * @param bankNames -A list of bankNames. See pwks
- * @param tolerance -A measure of the maxiumum distance a peak's h value,k
- *value, or l value is
- * from an integer to be considered indexed.
- * @return The associated workspace
- *
- * NOTE: This method could be used if this FitFunction is part of a composite
- *function, but an Xstart
- * and Xend for each composite is needed and may be difficult to
- *determine.
- */
- static DataObjects::Workspace2D_sptr
- calcWorkspace(DataObjects::PeaksWorkspace_sptr &pwks,
- std::vector<std::string> &bankNames, double tolerance);
-
- /**
- * Creates a new peak, matching the old peak except for a different
- *instrument.
- *
- * The Time of flightis the same except offset by T0. L0 should be the L0
- *for the new instrument.
- * It is added as a parameter in case the instrument will have the initial
- *flight path adjusted later.
- * NOTE: the wavelength is changed.
- *
- * @param peak_old - The old peak
- * @param instrNew -The new instrument
- * @param T0 :
- * @param L0 :
- * @return The new peak with the new instrument( adjusted with the
- *parameters) and time adjusted.
- */
- static DataObjects::Peak createNewPeak(const Geometry::IPeak &peak_old,
- Geometry::Instrument_sptr instrNew,
- double T0, double L0);
-
-protected:
- void init() override;
-
- /**
- * Checks for out of bounds values , peaksWorkspace status
- */
- void Check(DataObjects::PeaksWorkspace_sptr &pkwsp, const double *xValues,
- const size_t nData, size_t &StartX, size_t &EndX) const;
-
- /**
- * Gets the new instrument by applying parameters values to the old
- *instrument
- * @param peak A peak. Only used to get an old instrument from the 1st
- *peak.
- *
- * @return A new instrument with the parameters applied.
- */
- Geometry::Instrument_sptr getNewInstrument(const Geometry::IPeak &peak) const;
-
-private:
- /**
- * Even though constrains are used. Often very illogical parameters have to be
- * processed.
- * This checks for these conditions.
- */
- double checkForNonsenseParameters() const;
-
- boost::shared_ptr<const Geometry::IComponent> findBank(std::string bankName);
+ const IFunction::Attribute &value) override;
- std::map<std::string, boost::shared_ptr<const Geometry::IComponent>>
- bankDetMap;
+ /// Move detectors with parameters
+ void moveDetector(
+ double x, double y, double z, double rotx, double roty, double rotz,
+ std::string detname, API::Workspace_sptr inputW) const;
- /**
- * Get the peaks workspace that was specified.
- */
- void getPeaks() const;
+private:
+ /// Call the appropriate load function
+ void load(const std::string &fname);
- mutable boost::shared_ptr<DataObjects::PeaksWorkspace> m_peaks;
+ /// Load the points from a Workspace
+ void loadWorkspace(const std::string &wsName) const;
- double a, b, c, alpha, beta, gamma;
- int NGroups;
- bool RotateCenters, SampleOffsets;
- double SampleX, SampleY, SampleZ;
+ /// Load the points from a Workspace
+ void loadWorkspace(boost::shared_ptr<API::Workspace> ws) const;
- std::string PeakName; //< SCDPanelErrors{PeakName} is name in the Analysis
- // Data Service where the PeaksWorkspace is stored
+ /// Clear all data
+ void clear() const;
- bool a_set, b_set, c_set, alpha_set, beta_set, gamma_set, PeakName_set,
- BankNames_set, startX_set, endX_set, NGroups_set, sampleX_set,
- sampleY_set, sampleZ_set;
+ /// Evaluate the function for a list of arguments and given scaling factor
+ void eval(double xrotate, double yrotate, double zrotate, double xshift, double yshift, double zshift, double *out,
+ const double *xValues, const size_t nData) const;
- double tolerance;
+ /// Fill in the workspace and bank names
+ void setupData() const;
- /// The OrientedLattice created from the parameters
- boost::shared_ptr<Geometry::UnitCell> m_unitCell;
+ /// The default value for the workspace index
+ static const int defaultIndexValue;
- std::string BankNames;
+ /// Temporary workspace holder
+ mutable boost::shared_ptr<API::Workspace> m_workspace;
- int m_startX; ///<start index in xValues array in functionMW. -1 use all.
- int m_endX; ///<end index in xValues array in functionMW. -1 use all.
+ /// Stores bank
+ mutable std::string m_bank;
- std::vector<std::string> m_attrNames;
+ /// Flag of completing data setup
+ mutable bool m_setupFinished;
- // ki-kf for Inelastic convention; kf-ki for Crystallography convention
- std::string convention;
};
-}
-}
+
+} // namespace Crystal
+} // namespace Mantid
#endif /*MANTID_CRYSTAL_SCDPANELERRORS_H_*/
diff --git a/Framework/Crystal/src/PeakHKLErrors.cpp b/Framework/Crystal/src/PeakHKLErrors.cpp
index 5e2849dc58b..90a0daced7a 100644
--- a/Framework/Crystal/src/PeakHKLErrors.cpp
+++ b/Framework/Crystal/src/PeakHKLErrors.cpp
@@ -10,13 +10,13 @@
#include "MantidGeometry/Crystal/IPeak.h"
#include "MantidAPI/ParamFunction.h"
#include "MantidCrystal/PeakHKLErrors.h"
-#include "MantidCrystal/SCDPanelErrors.h"
-
+#include "MantidAPI/AnalysisDataService.h"
#include <boost/math/special_functions/round.hpp>
using namespace Mantid::DataObjects;
using namespace Mantid::API;
using namespace Mantid::Kernel;
+using namespace Mantid::Kernel::Units;
using Mantid::Geometry::CompAssembly;
using Mantid::Geometry::IObjComponent_const_sptr;
using Mantid::Geometry::IComponent_const_sptr;
@@ -378,7 +378,7 @@ void PeakHKLErrors::function1D(double *out, const double *xValues,
int runNum = peak_old.getRunNumber();
std::string runNumStr = std::to_string(runNum);
Peak peak =
- SCDPanelErrors::createNewPeak(peak_old, instNew, 0, peak_old.getL1());
+ createNewPeak(peak_old, instNew, 0, peak_old.getL1());
size_t N = OptRuns.find("/" + runNumStr + "/");
if (N < OptRuns.size()) {
@@ -471,7 +471,7 @@ void PeakHKLErrors::functionDeriv1D(Jacobian *out, const double *xValues,
int peakNum = boost::math::iround(xValues[i]);
IPeak &peak_old = Peaks->getPeak(peakNum);
Peak peak =
- SCDPanelErrors::createNewPeak(peak_old, instNew, 0, peak_old.getL1());
+ createNewPeak(peak_old, instNew, 0, peak_old.getL1());
int runNum = peak_old.getRunNumber();
std::string runNumStr = std::to_string(runNum);
@@ -623,5 +623,38 @@ void PeakHKLErrors::functionDeriv1D(Jacobian *out, const double *xValues,
}
}
}
+
+Peak PeakHKLErrors::createNewPeak(const Geometry::IPeak &peak_old,
+ Geometry::Instrument_sptr instrNew,
+ double T0, double L0) {
+ Geometry::Instrument_const_sptr inst = peak_old.getInstrument();
+ if (inst->getComponentID() != instrNew->getComponentID()) {
+ g_log.error("All peaks must have the same instrument");
+ throw std::invalid_argument("All peaks must have the same instrument");
+ }
+
+ double T = peak_old.getTOF() + T0;
+
+ int ID = peak_old.getDetectorID();
+
+ Kernel::V3D hkl = peak_old.getHKL();
+ // peak_old.setDetectorID(ID); //set det positions
+ Peak peak(instrNew, ID, peak_old.getWavelength(), hkl,
+ peak_old.getGoniometerMatrix());
+
+ Wavelength wl;
+
+ wl.initialize(L0, peak.getL2(), peak.getScattering(), 0,
+ peak_old.getInitialEnergy(), 0.0);
+
+ peak.setWavelength(wl.singleFromTOF(T));
+ peak.setIntensity(peak_old.getIntensity());
+ peak.setSigmaIntensity(peak_old.getSigmaIntensity());
+ peak.setRunNumber(peak_old.getRunNumber());
+ peak.setBinCount(peak_old.getBinCount());
+
+ //!!!peak.setDetectorID(ID);
+ return peak;
+}
}
}
diff --git a/Framework/Crystal/src/SCDCalibratePanels.cpp b/Framework/Crystal/src/SCDCalibratePanels.cpp
index 83b660caf27..e26b4df2994 100644
--- a/Framework/Crystal/src/SCDCalibratePanels.cpp
+++ b/Framework/Crystal/src/SCDCalibratePanels.cpp
@@ -9,10 +9,12 @@
#include "MantidAPI/IFunction.h"
#include "MantidAPI/FunctionFactory.h"
#include "MantidAPI/IFunction1D.h"
+#include "MantidCrystal/SCDPanelErrors.h"
#include <fstream>
#include "MantidGeometry/Crystal/IndexingUtils.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidGeometry/Crystal/ReducedCell.h"
+#include <boost/math/special_functions/round.hpp>
#include <Poco/File.h>
#include <sstream>
@@ -605,9 +607,13 @@ void SCDCalibratePanels::exec() {
}
}
}
-
- int nPeaks = peaksWs->getNumberPeaks();
- std::string cell_type = getProperty("CellType");
+ IAlgorithm_sptr ub_alg;
+ try {
+ ub_alg = createChildAlgorithm("CalculateUMatrix", -1, -1, false);
+ } catch (Exception::NotFoundError &) {
+ g_log.error("Can't locate CalculateUMatrix algorithm");
+ throw;
+ }
double a = getProperty("a");
double b = getProperty("b");
double c = getProperty("c");
@@ -624,656 +630,80 @@ void SCDCalibratePanels::exec() {
alpha = latt.alpha();
beta = latt.beta();
gamma = latt.gamma();
- } /*else {
- boost::shared_ptr<Algorithm> alg =
- createChildAlgorithm("FindUBUsingLatticeParameters", .2, .9, true);
- alg->setProperty("PeaksWorkspace", peaksWs);
- alg->setProperty("a", a);
- alg->setProperty("b", b);
- alg->setProperty("c", c);
- alg->setProperty("alpha", alpha);
- alg->setProperty("beta", beta);
- alg->setProperty("gamma", gamma);
- alg->setProperty("NumInitial", 15);
- alg->setProperty("Tolerance", 0.12);
- alg->executeAsChildAlg();
- }*/
- double tolerance = getProperty("tolerance");
-
- string DetCalFileName = getProperty("DetCalFilename");
- if (Poco::File(DetCalFileName).exists())
- Poco::File(DetCalFileName).remove();
-
- bool useL0 = getProperty("useL0");
- bool useTimeOffset = getProperty("useTimeOffset");
- int nIter = 1;
- if (useL0 || useTimeOffset)
- nIter = 2;
- bool use_PanelWidth = getProperty("usePanelWidth");
- bool use_PanelHeight = getProperty("usePanelHeight");
- bool use_PanelPosition = getProperty("usePanelPosition");
- bool use_PanelOrientation = getProperty("usePanelOrientation");
- double SampleXoffset = getProperty("SampleXoffset");
- double SampleYoffset = getProperty("SampleYoffset");
- double SampleZoffset = getProperty("SampleZoffset");
-
- string Grouping = getProperty("PanelGroups");
- string bankPrefix = getProperty("PanelNamePrefix");
- string bankingCode = getProperty("Grouping");
-
- //----------------- Set Up Bank Name Vectors -------------------------
- set<string, compareBanks> AllBankNames;
- for (int i = 0; i < nPeaks; ++i)
- AllBankNames.insert(peaksWs->getPeak(i).getBankName());
-
- vector<vector<string>> Groups;
- CalculateGroups(AllBankNames, Grouping, bankPrefix, bankingCode, Groups);
-
- //----------------- Calculate & Create Calculated vs Theoretical
- // workspaces------------------,);
- int nGroups = static_cast<int>(AllBankNames.size());
- MatrixWorkspace_sptr ColWksp =
- Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", nGroups,
- nPeaks, nPeaks);
- MatrixWorkspace_sptr RowWksp =
- Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", nGroups,
- nPeaks, nPeaks);
- MatrixWorkspace_sptr TofWksp =
- Mantid::API::WorkspaceFactory::Instance().create("Workspace2D", nGroups,
- nPeaks, nPeaks);
- nGroups = static_cast<int>(Groups.size());
-
- double chisqSum = 0;
- int NDofSum = 0;
- if (!GoodStart(peaksWs, a, b, c, alpha, beta, gamma, tolerance)) {
- g_log.warning() << "**** Indexing is NOT compatible with given lattice "
- "parameters******\n";
- g_log.warning()
- << " Index with conventional orientation matrix???\n";
}
+ ub_alg->setProperty("PeaksWorkspace", peaksWs);
+ ub_alg->setProperty("a", a);
+ ub_alg->setProperty("b", b);
+ ub_alg->setProperty("c", c);
+ ub_alg->setProperty("alpha", alpha);
+ ub_alg->setProperty("beta", beta);
+ ub_alg->setProperty("gamma", gamma);
+ ub_alg->executeAsChildAlg();
- //----------- Initialize peaksWorkspace, initial parameter values
- // etc.---------
- boost::shared_ptr<const Instrument> instrument =
- peaksWs->getPeak(0).getInstrument();
-
- // Time offset from property
- const API::Run &run = peaksWs->run();
- double T0 = 0.0;
- if (run.hasProperty("T0")) {
- Kernel::Property *prop = run.getProperty("T0");
- T0 = boost::lexical_cast<double, std::string>(prop->value());
- if (T0 != 0) {
- g_log.notice() << "T0 = " << T0 << '\n';
+ for (int i = int(peaksWs->getNumberPeaks()) - 1; i >= 0; --i) {
+ Peak peak = peaksWs->getPeak(i);
+ if (peak.getHKL() == V3D(0,0,0) || peak.getBankName() != "bank13") {
+ peaksWs->removePeak(i);
}
}
- if ((string)getProperty("PreProcessInstrument") ==
- "C)Apply a LoadParameter.xml type file")
- T0 = getProperty("InitialTimeOffset"); //!*****
-
- double L0 = peaksWs->getPeak(0).getL1();
- std::ofstream outfile(DetCalFileName);
- g_log.debug() << "Initial L0,T0=" << L0 << "," << T0 << '\n';
- std::vector<std::string> detcal;
-
- for (auto iIter = 0; iIter != nIter; ++iIter) {
- detcal.clear();
- PARALLEL_FOR1(peaksWs)
- for (int iGr = 0; iGr < nGroups; iGr++) {
- PARALLEL_START_INTERUPT_REGION
- auto group = Groups.begin() + iGr;
- vector<string> banksVec;
- for (auto &bankName : *group) {
- banksVec.push_back(bankName);
- }
- //------------------ Set Up Workspace for IFitFunction Fit---------------
- vector<int> bounds;
- Workspace2D_sptr ws = calcWorkspace(peaksWs, banksVec, tolerance, bounds);
- double T0_bank = T0;
- double L0_bank = L0;
- boost::shared_ptr<const Instrument> PreCalibinstrument =
- GetNewCalibInstrument(instrument,
- (string)getProperty("PreProcessInstrument"),
- (string)getProperty("PreProcFilename"), T0_bank,
- L0_bank, banksVec);
- V3D samplePos =
- peaksWs->getPeak(0).getInstrument()->getSample()->getPos();
-
- string PeakWSName = getPropertyValue("PeakWorkspace");
- if (PeakWSName.length() < 1) {
- PeakWSName = "xxx";
- AnalysisDataService::Instance().addOrReplace("xxx", peaksWs);
- }
-
- int NGroups = 1; //(int)Groups.size();
- double detWidthScale0, detHeightScale0, Xoffset0, Yoffset0, Zoffset0,
- Xrot0, Yrot0, Zrot0;
-
- //------------------- For each Group set up Function,
- //--------------------------
- //---------------Ties, and Constraint Properties for Fit
- // algorithm--------------------
-
- // set up the string for specifying groups
- string BankNameString;
- // for (auto group = Groups.begin(); group != Groups.end(); ++group) {
- // if (group != Groups.begin())
- // BankNameString += "!";
- for (auto bank = group->begin(); bank != group->end(); ++bank) {
- if (bank != group->begin())
- BankNameString += "/";
-
- BankNameString += (*bank);
- }
- //}
-
- int RotGroups = 0;
- if (getProperty("RotateCenters"))
- RotGroups = 1;
- int SampOffsets = 0;
- if (getProperty("AllowSampleShift"))
- SampOffsets = 1;
-
- // first round of function setup
- IFunction_sptr iFunc =
- FunctionFactory::Instance().createFunction("SCDPanelErrors");
- iFunc->setAttributeValue("PeakWorkspaceName", PeakWSName);
- iFunc->setAttributeValue("a", a);
- iFunc->setAttributeValue("b", b);
- iFunc->setAttributeValue("c", c);
- iFunc->setAttributeValue("alpha", alpha);
- iFunc->setAttributeValue("beta", beta);
- iFunc->setAttributeValue("gamma", gamma);
- iFunc->setAttributeValue("NGroups", NGroups);
- iFunc->setAttributeValue("BankNames", BankNameString);
- iFunc->setAttributeValue("startX", -1);
- iFunc->setAttributeValue("endX", -1);
- iFunc->setAttributeValue("RotateCenters", RotGroups);
- iFunc->setAttributeValue("SampleOffsets", SampOffsets);
- iFunc->setParameter("l0", L0_bank);
- iFunc->setParameter("t0", T0_bank);
-
- double maxXYOffset = getProperty("MaxPositionChange_meters");
-
- boost::shared_ptr<const RectangularDetector> bank_rect;
-
- string name = group->front();
- boost::shared_ptr<const IComponent> bank_cmp =
- instrument->getComponentByName(name);
- bank_rect =
- boost::dynamic_pointer_cast<const RectangularDetector>(bank_cmp);
-
- if (!bank_rect) {
- g_log.error("No Rectangular detector bank " + banksVec[0] +
- " in instrument");
- throw invalid_argument("No Rectangular detector bank " + banksVec[0] +
- " in instrument");
- }
-
- // if( it1 == (*itv).begin())
- CalcInitParams(bank_rect, instrument, PreCalibinstrument, detWidthScale0,
- detHeightScale0, Xoffset0, Yoffset0, Zoffset0, Xrot0,
- Yrot0, Zrot0);
-
- // --- set Function property ----------------------
- iFunc->setParameter("f0_detWidthScale", detWidthScale0);
- iFunc->setParameter("f0_detHeightScale", detHeightScale0);
- iFunc->setParameter("f0_Xoffset", Xoffset0);
- iFunc->setParameter("f0_Yoffset", Yoffset0);
- iFunc->setParameter("f0_Zoffset", Zoffset0);
- iFunc->setParameter("f0_Xrot", Xrot0);
- iFunc->setParameter("f0_Yrot", Yrot0);
- iFunc->setParameter("f0_Zrot", Zrot0);
-
- set<string> MyBankNames;
- for (int i = 0; i < nPeaks; ++i)
- MyBankNames.insert(banksVec[0]);
- int startX = bounds[0];
- int endXp1 = bounds[group->size()];
- if (endXp1 - startX < 13) {
- g_log.error() << "Bank Group " << BankNameString
- << " does not have enough peaks for fitting\n";
- saveIsawDetCal(instrument, MyBankNames, T0_bank,
- DetCalFileName + std::to_string(iGr));
- continue;
- }
-
- //---------- setup ties ----------------------------------
- tie(iFunc, !use_PanelWidth, "f0_detWidthScale", detWidthScale0);
- tie(iFunc, !use_PanelHeight, "f0_detHeightScale", detHeightScale0);
- tie(iFunc, !use_PanelPosition, "f0_Xoffset", Xoffset0);
- tie(iFunc, !use_PanelPosition, "f0_Yoffset", Yoffset0);
- tie(iFunc, !use_PanelPosition, "f0_Zoffset", Zoffset0);
- tie(iFunc, !use_PanelOrientation, "f0_Xrot", Xrot0);
- tie(iFunc, !use_PanelOrientation, "f0_Yrot", Yrot0);
- tie(iFunc, !use_PanelOrientation, "f0_Zrot", Zrot0);
-
- //--------------- setup constraints ------------------------------
- constrain(iFunc, "l0", (MIN_DET_HW_SCALE * L0_bank),
- (MAX_DET_HW_SCALE * L0_bank));
- constrain(iFunc, "t0", -5., 5.);
-
- constrain(iFunc, "f0_detWidthScale", MIN_DET_HW_SCALE * detWidthScale0,
- MAX_DET_HW_SCALE * detWidthScale0);
- constrain(iFunc, "f0_detHeightScale", MIN_DET_HW_SCALE * detHeightScale0,
- MAX_DET_HW_SCALE * detHeightScale0);
- constrain(iFunc, "f0_Xoffset", -1. * maxXYOffset + Xoffset0,
- maxXYOffset + Xoffset0);
- constrain(iFunc, "f0_Yoffset", -1. * maxXYOffset + Yoffset0,
- maxXYOffset + Yoffset0);
- constrain(iFunc, "f0_Zoffset", -1. * maxXYOffset + Zoffset0,
- maxXYOffset + Zoffset0);
-
- double MaxRotOffset = getProperty("MaxRotationChangeDegrees");
- constrain(iFunc, "f0_Xrot", -1. * MaxRotOffset, MaxRotOffset);
- constrain(iFunc, "f0_Yrot", -1. * MaxRotOffset, MaxRotOffset);
- constrain(iFunc, "f0_Zrot", -1. * MaxRotOffset, MaxRotOffset);
- //} // for vector< string > in Groups
-
- // Function supports setting the sample position even when it isn't be
- // refined
- iFunc->setAttributeValue("SampleX", samplePos.X() + SampleXoffset);
- iFunc->setAttributeValue("SampleY", samplePos.Y() + SampleYoffset);
- iFunc->setAttributeValue("SampleZ", samplePos.Z() + SampleZoffset);
-
- // Constraints for sample offsets
- if (getProperty("AllowSampleShift")) {
- maxXYOffset = getProperty("MaxSamplePositionChangeMeters");
- constrain(iFunc, "SampleX", samplePos.X() + SampleXoffset - maxXYOffset,
- samplePos.X() + SampleXoffset + maxXYOffset);
- constrain(iFunc, "SampleY", samplePos.Y() + SampleYoffset - maxXYOffset,
- samplePos.Y() + SampleYoffset + maxXYOffset);
- constrain(iFunc, "SampleZ", samplePos.Z() + SampleZoffset - maxXYOffset,
- samplePos.Z() + SampleZoffset + maxXYOffset);
- } else {
- tie(iFunc, true, "SampleX", samplePos.X() + SampleXoffset);
- tie(iFunc, true, "SampleY", samplePos.Y() + SampleYoffset);
- tie(iFunc, true, "SampleZ", samplePos.Z() + SampleZoffset);
- }
-
- tie(iFunc, !useL0, "l0", L0_bank);
- tie(iFunc, !useTimeOffset, "t0", T0_bank);
-
- //--------------------- Set up Fit Algorithm and
- // Execute-------------------
- boost::shared_ptr<Algorithm> fit_alg =
- createChildAlgorithm("Fit", .2, .9, true);
-
- if (!fit_alg)
- throw invalid_argument("Cannot find Fit algorithm");
- fit_alg->initialize();
-
- int Niterations = getProperty("NumIterations");
- std::string minimizerError = getProperty("MinimizerError");
- std::string minimizer = getProperty("Minimizer");
- fit_alg->setProperty("Function", iFunc);
- fit_alg->setProperty("MaxIterations", Niterations);
- fit_alg->setProperty("InputWorkspace", ws);
- fit_alg->setProperty("Output", "out");
- fit_alg->setProperty("CreateOutput", true);
- fit_alg->setProperty("CalcErrors", false);
- fit_alg->setPropertyValue("Minimizer", minimizer + ",AbsError=" +
- minimizerError + ",RelError=" +
- minimizerError);
- fit_alg->executeAsChildAlg();
-
- PARALLEL_CRITICAL(afterFit) {
- MatrixWorkspace_sptr fitWS = fit_alg->getProperty("OutputWorkspace");
- // AnalysisDataService::Instance().addOrReplace("out"+boost::lexical_cast<std::string>(iGr),
- // fitWS);
- g_log.debug() << "Finished executing Fit algorithm\n";
- string OutputStatus = fit_alg->getProperty("OutputStatus");
- g_log.notice() << BankNameString << " Output Status=" << OutputStatus
- << "\n";
-
- //--------------------- Get and Process Results -----------------------
- double chisq = fit_alg->getProperty("OutputChi2overDoF");
- chisqSum += chisq;
-
- if (chisq > 1) {
- g_log.warning() << "************* This is a large chi squared value "
- "************\n";
- g_log.warning()
- << " the indexing may have been using an incorrect\n";
- g_log.warning() << " orientation matrix, instrument geometry or "
- "goniometer info\n";
- }
- ITableWorkspace_sptr RRes = fit_alg->getProperty("OutputParameters");
- vector<double> params;
- vector<double> errs;
- vector<string> names;
- double sigma = sqrt(chisq);
-
- if (chisq < 0 || chisq != chisq)
- sigma = -1;
- string fieldBaseNames =
- ";l0;t0;detWidthScale;detHeightScale;Xoffset;Yoffset;"
- "Zoffset;Xrot;Yrot;Zrot;";
- if (getProperty("AllowSampleShift"))
- fieldBaseNames += "SampleX;SampleY;SampleZ;";
- for (size_t prm = 0; prm < RRes->rowCount(); ++prm) {
- string namee = RRes->getRef<string>("Name", prm);
- size_t dotPos = namee.find('_');
- if (dotPos >= namee.size())
- dotPos = 0;
- else
- dotPos++;
- string Field = namee.substr(dotPos);
- size_t FieldNum = fieldBaseNames.find(";" + Field + ";");
- if (FieldNum > fieldBaseNames.size())
- continue;
- if (dotPos != 0) {
- int col = atoi(namee.substr(1, dotPos).c_str());
- if (col < 0 || col >= NGroups)
- continue;
- }
- names.push_back(namee);
- params.push_back(RRes->getRef<double>("Value", prm));
- double err = RRes->getRef<double>("Error", prm);
- errs.push_back(sigma * err);
- }
-
- //------------------- Report chi^2 value --------------------
- int nVars = 8; // NGroups;
-
- if (!use_PanelWidth)
- nVars--;
- if (!use_PanelHeight)
- nVars--;
- if (!use_PanelPosition)
- nVars -= 3;
- if (!use_PanelOrientation)
- nVars -= 3;
- nVars *= NGroups;
- nVars += 2;
-
- if (!useL0)
- nVars--;
- if (!useTimeOffset)
- nVars--;
-
- // g_log.notice() << " nVars=" <<nVars<< '\n';
- int NDof = (static_cast<int>(ws->dataX(0).size()) - nVars);
- NDofSum = +NDof;
-
- map<string, double> result;
-
- for (size_t i = 0; i < min<size_t>(params.size(), names.size()); ++i) {
- result[names[i]] = params[i];
- }
-
- g_log.notice() << BankNameString << " ChiSqoverDoF =" << chisq
- << " NDof =" << NDof << " l0 = " << result["l0"]
- << " T0 = " << result["t0"]
- << " peaks = " << endXp1 - startX << "\n";
-
- //--------------------- Create Result Table Workspace-------------------
- this->progress(.92, "Creating Results table");
- createResultWorkspace(nGroups, iGr + 1, names, params, errs);
-
- //---------------- Create new instrument with ------------------------
- //--------------new parameters to SAVE to files---------------------
-
- auto pmap = boost::make_shared<ParameterMap>();
- boost::shared_ptr<const ParameterMap> pmapOld =
- instrument->getParameterMap();
- boost::shared_ptr<const Instrument> NewInstrument =
- boost::make_shared<Instrument>(instrument->baseInstrument(), pmap);
-
- boost::shared_ptr<const RectangularDetector> bank_rect;
- double rotx, roty, rotz;
-
- rotx = result["f0_Xrot"];
- roty = result["f0_Yrot"];
- rotz = result["f0_Zrot"];
-
- Quat newRelRot = Quat(rotx, V3D(1, 0, 0)) * Quat(roty, V3D(0, 1, 0)) *
- Quat(rotz, V3D(0, 0, 1)); //*RelRot;
-
- FixUpBankParameterMap((banksVec), NewInstrument,
- V3D(result["f0_Xoffset"], result["f0_Yoffset"],
- result["f0_Zoffset"]),
- newRelRot, result["f0_detWidthScale"],
- result["f0_detHeightScale"], pmapOld,
- getProperty("RotateCenters"));
-
- //} // For @ group
-
- V3D sampPos(
- NewInstrument->getSample()->getPos()); // should be (0,0,0)???
- if (getProperty("AllowSampleShift"))
- sampPos =
- V3D(result["SampleX"], result["SampleY"], result["SampleZ"]);
+ int nPeaks = peaksWs->getNumberPeaks();
- FixUpSourceParameterMap(NewInstrument, result["l0"], sampPos, pmapOld);
+ MatrixWorkspace_sptr q3DWS = boost::dynamic_pointer_cast<MatrixWorkspace>(
+ API::WorkspaceFactory::Instance().create("Workspace2D", 1,
+ nPeaks*3, nPeaks*3));
- //---------------------- Save new instrument to DetCal--------------
- this->progress(.94, "Saving detcal file");
- saveIsawDetCal(NewInstrument, MyBankNames, result["t0"],
- DetCalFileName + std::to_string(iGr));
- }
- PARALLEL_END_INTERUPT_REGION
- }
- PARALLEL_CHECK_INTERUPT_REGION
- std::vector<double> l0vec, t0vec;
- std::string line, seven;
- for (int iGr = 0; iGr < nGroups; iGr++) {
- std::ifstream infile(DetCalFileName + std::to_string(iGr));
- while (std::getline(infile, line)) {
- if (iGr == 0 && iIter == nIter - 1) {
- if (line[0] == '#' || line[0] == '6')
- outfile << line << "\n";
- }
- if (line[0] == '7') {
- double L0bank, T0bank;
- std::stringstream(line) >> seven >> L0bank >> T0bank;
- l0vec.push_back(L0bank);
- t0vec.push_back(T0bank);
- } else if (line[0] == '5')
- detcal.push_back(line);
- }
- infile.close();
- if (Poco::File(DetCalFileName + std::to_string(iGr)).exists())
- Poco::File(DetCalFileName + std::to_string(iGr)).remove();
- }
- if (useL0) {
- removeOutliers(l0vec);
- Statistics stats = getStatistics(l0vec);
- L0 = stats.mean * 0.01; // cm when read from file
- useL0 = false;
- }
- if (useTimeOffset) {
- removeOutliers(t0vec);
- Statistics stats = getStatistics(t0vec);
- T0 = stats.mean;
- // set T0 in the run parameters
- API::Run &m_run = peaksWs->mutableRun();
- m_run.addProperty<double>("T0", T0, true);
- useTimeOffset = false;
- }
- }
- L0 *= 100; // ISAW uses cm
- outfile << "7 " << std::setprecision(4) << std::fixed << L0;
- outfile << std::setw(13) << std::setprecision(3) << T0 << '\n';
- outfile << "4 DETNUM NROWS NCOLS WIDTH HEIGHT DEPTH DETD CenterX "
- " CenterY CenterZ BaseX BaseY BaseZ UpX UpY "
- " UpZ\n";
- for (vector<std::string>::const_iterator itdet = detcal.begin();
- itdet != detcal.end(); ++itdet)
- outfile << *itdet << "\n";
- outfile.close();
-
- setProperty("ColWorkspace", ColWksp);
- setProperty("RowWorkspace", RowWksp);
- setProperty("TofWorkspace", TofWksp);
- setProperty("ChiSqOverDOF", chisqSum);
- setProperty("DOF", NDofSum);
-
- set<string> bankNames;
- for (int i = 0; i < nPeaks; ++i) {
- IPeak &peak = peaksWs->getPeak(i);
- instrument = peak.getInstrument();
- LoadISawDetCal(instrument, bankNames, T0, L0, DetCalFileName, "bank");
- peak.setInstrument(instrument);
- peak.setDetectorID(peak.getDetectorID());
- }
+ auto &outSpec = q3DWS->getSpectrum(0);
+ MantidVec &yVec = outSpec.dataY();
+ MantidVec &xVec = outSpec.dataX();
- //-----------------------Save new instrument to xml(for LoadParameterFile)
- // files----------
- this->progress(.96, "Saving xml param file");
- string XmlFileName = getProperty("XmlFilename");
- saveXmlFile(XmlFileName, Groups, peaksWs->getInstrument());
-
- IFunction_sptr fn =
- FunctionFactory::Instance().createFunction("LatticeFunction");
- fn->setAttributeValue("LatticeSystem", cell_type);
- fn->addTies("ZeroShift=0.0");
- fn->setParameter("a", a);
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::MONOCLINIC() ||
- cell_type == ReducedCell::ORTHORHOMBIC()) {
- fn->setParameter("b", b);
- }
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::TETRAGONAL() ||
- cell_type == ReducedCell::ORTHORHOMBIC() ||
- cell_type == ReducedCell::HEXAGONAL() ||
- cell_type == ReducedCell::MONOCLINIC()) {
- fn->setParameter("c", c);
- }
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::RHOMBOHEDRAL()) {
- fn->setParameter("Alpha", alpha);
- }
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::MONOCLINIC()) {
- fn->setParameter("Beta", beta);
- }
- if (cell_type == ReducedCell::TRICLINIC()) {
- fn->setParameter("Gamma", gamma);
+ OrientedLattice lattice = peaksWs->mutableSample().getOrientedLattice();
+ for (int i = 0; i < nPeaks; i++) {
+ Peak peak = peaksWs->getPeak(i);
+ V3D hkl = V3D(boost::math::iround(peak.getH()), boost::math::iround(peak.getK()),
+ boost::math::iround(peak.getL()));
+ V3D Q2 = lattice.qFromHKL(hkl);
+ xVec[i*3] = i*3;
+ xVec[i*3+1] = i*3+1;
+ xVec[i*3+2] = i*3+2;
+ yVec[i*3] = Q2.X();
+ yVec[i*3+1] = Q2.Y();
+ yVec[i*3+2] = Q2.Z();
}
- IAlgorithm_sptr fit =
- Mantid::API::AlgorithmFactory::Instance().create("Fit", -1);
- fit->initialize();
- fit->setChild(true);
- fit->setLogging(false);
- fit->setProperty("Function", fn);
- fit->setProperty("InputWorkspace", peaksWs);
- fit->setProperty("CostFunction", "Unweighted least squares");
- fit->setProperty("CreateOutput", true);
- fit->execute();
-
- IAlgorithm_sptr ub_alg;
+ IAlgorithm_sptr fit_alg;
try {
- ub_alg = createChildAlgorithm("CalculateUMatrix", -1, -1, false);
+ fit_alg = createChildAlgorithm("Fit", -1, -1, false);
} catch (Exception::NotFoundError &) {
- g_log.error("Can't locate CalculateUMatrix algorithm");
+ g_log.error("Can't locate Fit algorithm");
throw;
}
+ std::ostringstream fun_str;
+ fun_str << "name=SCDPanelErrors,Workspace="<<peaksWs->name()<<",Bank=bank13";
+ fit_alg->setPropertyValue("Function", fun_str.str());
+ fit_alg->setProperty("InputWorkspace", q3DWS);
+ fit_alg->setProperty("CreateOutput", true);
+ fit_alg->setProperty("Output", "fit");
+ fit_alg->setProperty("MaxIterations", 0);
+ fit_alg->executeAsChildAlg();
+ MatrixWorkspace_sptr fitWS = fit_alg->getProperty("OutputWorkspace");
+ AnalysisDataService::Instance().addOrReplace("fit", fitWS);
+ ITableWorkspace_sptr paramsWS = fit_alg->getProperty("OutputParameters");
+ AnalysisDataService::Instance().addOrReplace("params", paramsWS);
+ double xShift = paramsWS->getRef<double>("Value",0);
+ double yShift = paramsWS->getRef<double>("Value",1);
+ double zShift = paramsWS->getRef<double>("Value",2);
+ double xRotate = paramsWS->getRef<double>("Value",3);
+ double yRotate = paramsWS->getRef<double>("Value",4);
+ double zRotate = paramsWS->getRef<double>("Value",5);
+ SCDPanelErrors det;
+ det.moveDetector(xShift, yShift, zShift, xRotate, yRotate, zRotate, "bank13", peaksWs);
- ub_alg->setProperty("PeaksWorkspace", peaksWs);
- ub_alg->setProperty("a", fn->getParameter("a"));
- double sigabc[6];
- sigabc[0] = fn->getError(0);
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::MONOCLINIC() ||
- cell_type == ReducedCell::ORTHORHOMBIC()) {
- ub_alg->setProperty("b", fn->getParameter("b"));
- sigabc[1] = fn->getError(1);
- } else {
- ub_alg->setProperty("b", fn->getParameter("a"));
- sigabc[1] = fn->getError(0);
- }
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::TETRAGONAL() ||
- cell_type == ReducedCell::ORTHORHOMBIC() ||
- cell_type == ReducedCell::HEXAGONAL() ||
- cell_type == ReducedCell::MONOCLINIC()) {
- ub_alg->setProperty("c", fn->getParameter("c"));
- sigabc[2] = fn->getError(2);
- } else {
- ub_alg->setProperty("c", fn->getParameter("a"));
- sigabc[2] = fn->getError(0);
- }
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::RHOMBOHEDRAL()) {
- ub_alg->setProperty("alpha", fn->getParameter("Alpha"));
- sigabc[3] = fn->getError(3);
- } else {
- ub_alg->setProperty("alpha", 90.0);
- sigabc[3] = 0.0;
- }
- if (cell_type == ReducedCell::TRICLINIC() ||
- cell_type == ReducedCell::MONOCLINIC()) {
- ub_alg->setProperty("beta", fn->getParameter("Beta"));
- sigabc[4] = fn->getError(4);
- } else {
- ub_alg->setProperty("beta", 90.0);
- sigabc[4] = 0.0;
- }
- if (cell_type == ReducedCell::TRICLINIC()) {
- ub_alg->setProperty("gamma", fn->getParameter("Gamma"));
- sigabc[5] = fn->getError(5);
- } else if (cell_type == ReducedCell::HEXAGONAL()) {
- ub_alg->setProperty("gamma", 120.0);
- sigabc[5] = 0.0;
- } else {
- ub_alg->setProperty("gamma", 90.0);
- sigabc[5] = 0.0;
- }
- ub_alg->executeAsChildAlg();
- OrientedLattice lattice = peaksWs->mutableSample().getOrientedLattice();
- DblMatrix UB = lattice.getUB();
- lattice.setError(sigabc[0], sigabc[1], sigabc[2], sigabc[3], sigabc[4],
- sigabc[5]);
- peaksWs->mutableSample().setOrientedLattice(&lattice);
- g_log.notice() << "Lattice after optimization: " << lattice << "\n";
-
- // We must sort the peaks
- // std::vector<std::pair<std::string, bool>> criteria;
- // criteria.push_back(std::pair<std::string, bool>("BankName", true));
- criteria.push_back(std::pair<std::string, bool>("RunNumber", true));
- criteria.push_back(std::pair<std::string, bool>("h", true));
- criteria.push_back(std::pair<std::string, bool>("k", true));
- criteria.push_back(std::pair<std::string, bool>("l", true));
- peaksWs->sort(criteria);
-
- // create table of theoretical vs calculated
- int bankLast = -1;
- int iSpectrum = -1;
- int icount = 0;
-
- for (int j = 0; j < nPeaks; ++j) {
- const Geometry::IPeak &peak = peaksWs->getPeak(j);
- string bankName = peak.getBankName();
- size_t k = bankName.find_last_not_of("0123456789");
- int bank = 0;
- if (k < bankName.length())
- bank = boost::lexical_cast<int>(bankName.substr(k + 1));
- if (bank != bankLast) {
- iSpectrum++;
- ColWksp->getSpectrum(iSpectrum).setSpectrumNo(specnum_t(bank));
- RowWksp->getSpectrum(iSpectrum).setSpectrumNo(specnum_t(bank));
- TofWksp->getSpectrum(iSpectrum).setSpectrumNo(specnum_t(bank));
- bankLast = bank;
- icount = 0;
- }
-
- try {
- V3D q_lab = (peak.getGoniometerMatrix() * UB) * peak.getHKL() * M_2_PI;
- Peak theoretical(peak.getInstrument(), q_lab);
- ColWksp->dataX(iSpectrum)[icount] = peak.getCol();
- ColWksp->dataY(iSpectrum)[icount] = theoretical.getCol();
- RowWksp->dataX(iSpectrum)[icount] = peak.getRow();
- RowWksp->dataY(iSpectrum)[icount] = theoretical.getRow();
- TofWksp->dataX(iSpectrum)[icount] = peak.getTOF();
- TofWksp->dataY(iSpectrum)[icount] = theoretical.getTOF();
- icount++;
- } catch (...) {
- // g_log.debug() << "Problem only in printing peaks\n";
- }
- }
+ string DetCalFileName = getProperty("DetCalFilename");
+ set<string> MyBankNames;
+ MyBankNames.insert("bank13");
+ Instrument_sptr inst = boost::const_pointer_cast<Instrument>(peaksWs->getInstrument());
+ saveIsawDetCal(inst, MyBankNames, 0.0,
+ DetCalFileName);
}
/**
@@ -1487,7 +917,7 @@ void SCDCalibratePanels::createResultWorkspace(const int numGroups,
Result->cell<double>(rowNum + 10 + nn, colNum) = errs[p];
}
- setProperty("ResultWorkspace", Result);
+ //setProperty("ResultWorkspace", Result);
}
/**
@@ -1501,7 +931,7 @@ void SCDCalibratePanels::createResultWorkspace(const int numGroups,
* @param filename -The name of the DetCal file to save the results to
*/
void SCDCalibratePanels::saveIsawDetCal(
- boost::shared_ptr<const Instrument> &instrument, set<string> &AllBankName,
+ boost::shared_ptr<Instrument> &instrument, set<string> &AllBankName,
double T0, string filename) {
// having a filename triggers doing the work
if (filename.empty())
@@ -1629,7 +1059,7 @@ void SCDCalibratePanels::init() {
"Path to an Mantid .xml description(for LoadParameterFile) file to "
"save.");
- declareProperty(
+ /*declareProperty(
Kernel::make_unique<WorkspaceProperty<ITableWorkspace>>(
"ResultWorkspace", "ResultWorkspace", Kernel::Direction::Output),
"Workspace of Results");
@@ -1647,15 +1077,15 @@ void SCDCalibratePanels::init() {
declareProperty(
Kernel::make_unique<WorkspaceProperty<MatrixWorkspace>>(
"TofWorkspace", "TofWorkspace", Kernel::Direction::Output),
- "Workspace comparing calculated and theoretical TOF of each peak.");
+ "Workspace comparing calculated and theoretical TOF of each peak.");*/
const string OUTPUTS("Outputs");
setPropertyGroup("DetCalFilename", OUTPUTS);
setPropertyGroup("XmlFilename", OUTPUTS);
- setPropertyGroup("ResultWorkspace", OUTPUTS);
+ /*setPropertyGroup("ResultWorkspace", OUTPUTS);
setPropertyGroup("ColWorkspace", OUTPUTS);
setPropertyGroup("RowWorkspace", OUTPUTS);
- setPropertyGroup("TofWorkspace", OUTPUTS);
+ setPropertyGroup("TofWorkspace", OUTPUTS);*/
//------------------------------------ Tolerance
// settings-------------------------
diff --git a/Framework/Crystal/src/SCDPanelErrors.cpp b/Framework/Crystal/src/SCDPanelErrors.cpp
index b1a2a2a48b6..47cc7b06634 100644
--- a/Framework/Crystal/src/SCDPanelErrors.cpp
+++ b/Framework/Crystal/src/SCDPanelErrors.cpp
@@ -1,711 +1,286 @@
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
#include "MantidCrystal/SCDPanelErrors.h"
-#include "MantidCrystal/SCDCalibratePanels.h"
+#include "MantidKernel/FileValidator.h"
+#include "MantidAPI/Algorithm.h"
+#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/FunctionFactory.h"
-#include "MantidAPI/WorkspaceFactory.h"
-#include "MantidGeometry/Crystal/OrientedLattice.h"
-#include "MantidGeometry/Crystal/IndexingUtils.h"
-#include "MantidKernel/ConfigService.h"
-
-using namespace Mantid::API;
-using namespace Mantid::DataObjects;
-using namespace Mantid::Geometry;
-using namespace Mantid::Kernel;
-using namespace Mantid::Kernel::Units;
-using std::invalid_argument;
-using std::logic_error;
-using std::map;
-using std::runtime_error;
-using std::setw;
-using std::string;
-using std::vector;
-
-// TODO: add const's and & to arguments. PeaksWorkspace_ptr should be const_ptr
+#include "MantidGeometry/Instrument.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
+#include "MantidDataObjects/Peak.h"
+#include <algorithm>
+#include <cmath>
+#include <fstream>
+#include <sstream>
namespace Mantid {
namespace Crystal {
-namespace {
-Kernel::Logger g_log("SCDPanelErrors");
-}
-
-DECLARE_FUNCTION(SCDPanelErrors)
-
-// Assumes UB from optimize UB maps hkl to qxyz/2PI. So conversion factors to an
-// from
-// UB ified q's are below.
-
-namespace { // anonymous namespace
-// static const double ONE_OVER_TWO_PI = 1. / M_2_PI;
-const string LATTICE_A("a");
-const string LATTICE_B("b");
-const string LATTICE_C("c");
-const string LATTICE_ALPHA("alpha");
-const string LATTICE_BETA("beta");
-const string LATTICE_GAMMA("gamma");
-const string BANK_NAMES("BankNames");
-const string NUM_GROUPS("NGroups");
-const string X_START("startX");
-const string X_END("endX");
-const string PEAKS_WKSP("PeakWorkspaceName");
-const string ROTATE_CEN("RotateCenters");
-const string SAMPLE_OFF("SampleOffsets");
-const string SAMPLE_X("SampleX");
-const string SAMPLE_Y("SampleY");
-const string SAMPLE_Z("SampleZ");
-}
-
-void initializeAttributeList(vector<string> &attrs) {
- attrs.clear();
- attrs.push_back(LATTICE_A);
- attrs.push_back(LATTICE_B);
- attrs.push_back(LATTICE_C);
- attrs.push_back(LATTICE_ALPHA);
- attrs.push_back(LATTICE_BETA);
- attrs.push_back(LATTICE_GAMMA);
- attrs.push_back(PEAKS_WKSP);
- attrs.push_back(BANK_NAMES);
- attrs.push_back(X_START);
- attrs.push_back(X_END);
- attrs.push_back(NUM_GROUPS);
- attrs.push_back(ROTATE_CEN);
- attrs.push_back(SAMPLE_OFF);
- attrs.push_back(SAMPLE_X);
- attrs.push_back(SAMPLE_Y);
- attrs.push_back(SAMPLE_Z);
-}
-
-SCDPanelErrors::SCDPanelErrors()
- : API::ParamFunction(), IFunction1D(), tolerance(.6), m_startX(-1),
- m_endX(-1) {
- convention = Kernel::ConfigService::Instance().getString("Q.convention");
- initializeAttributeList(m_attrNames);
-
- SampleX = SampleY = SampleZ = 0.;
-
- a_set = b_set = c_set = alpha_set = beta_set = gamma_set = PeakName_set =
- BankNames_set = endX_set = startX_set = NGroups_set = sampleX_set =
- sampleY_set = sampleZ_set = false;
-
- // g_log.setLevel(7);
-
- BankNames = "";
-
- PeakName = "";
-
- a = b = c = alpha = beta = gamma = 0;
-
- NGroups = 1;
- RotateCenters = false;
- SampleOffsets = false;
-}
-
-size_t SCDPanelErrors::nAttributes() const { return m_attrNames.size(); }
-
-std::vector<std::string> SCDPanelErrors::getAttributeNames() const {
- return m_attrNames;
-}
-
-IFunction::Attribute
-SCDPanelErrors::getAttribute(const std::string &attName) const {
- if (!hasAttribute(attName))
- throw std::invalid_argument("Not a valid attribute name \"" + attName +
- "\"");
- if (attName == LATTICE_A)
- return Attribute(a);
- if (attName == LATTICE_B)
- return Attribute(b);
- if (attName == LATTICE_C)
- return Attribute(c);
- if (attName == LATTICE_ALPHA)
- return Attribute(alpha);
- if (attName == LATTICE_BETA)
- return Attribute(beta);
- if (attName == LATTICE_GAMMA)
- return Attribute(gamma);
- if (attName == (BANK_NAMES))
- return Attribute(BankNames);
- else if (attName == PEAKS_WKSP)
- return Attribute(PeakName);
- else if (attName == NUM_GROUPS)
- return Attribute(NGroups);
- else if (attName == ROTATE_CEN) {
- if (RotateCenters)
- return Attribute(1);
- else
- return Attribute(0);
- } else if (attName == SAMPLE_OFF) {
- if (SampleOffsets)
+using namespace CurveFitting;
- return Attribute(1);
- else
- return Attribute(0);
- } else if (attName == X_START)
- return Attribute(static_cast<int>(m_startX));
- else if (attName == X_END)
- return Attribute(static_cast<int>(m_endX));
- else if (attName == SAMPLE_X)
- return Attribute(SampleX);
- else if (attName == SAMPLE_Y)
- return Attribute(SampleY);
- else if (attName == SAMPLE_Z)
- return Attribute(SampleZ);
+using namespace Kernel;
- throw std::invalid_argument("Not a valid attribute name \"" + attName + "\"");
-}
+using namespace API;
-bool SCDPanelErrors::hasAttribute(const std::string &attName) const {
- return (std::find(m_attrNames.begin(), m_attrNames.end(), attName) !=
- m_attrNames.end());
+namespace {
+/// static logger
+Logger g_log("SCDPanelErrors");
}
-SCDPanelErrors::SCDPanelErrors(DataObjects::PeaksWorkspace_sptr &pwk,
- std::string &Component_name, double ax,
- double bx, double cx, double alphax,
- double betax, double gammax, double tolerance1)
- : API::ParamFunction(), IFunction1D() {
- initializeAttributeList(m_attrNames);
-
- m_peaks = pwk;
- BankNames = Component_name;
-
- tolerance = tolerance1;
- NGroups = 1;
- a_set = b_set = c_set = alpha_set = beta_set = gamma_set = PeakName_set =
- BankNames_set = endX_set = startX_set = NGroups_set = false;
-
- setAttribute(LATTICE_A, Attribute(ax));
- setAttribute(LATTICE_B, Attribute(bx));
- setAttribute(LATTICE_C, Attribute(cx));
- setAttribute(LATTICE_ALPHA, Attribute(alphax));
- setAttribute(LATTICE_BETA, Attribute(betax));
- setAttribute(LATTICE_GAMMA, Attribute(gammax));
+DECLARE_FUNCTION(SCDPanelErrors)
- setAttribute(PEAKS_WKSP, Attribute("xxx"));
- setAttribute(BANK_NAMES, Attribute(Component_name));
- setAttribute(X_START, Attribute(-1));
- setAttribute(X_END, Attribute(-1));
- setAttribute(ROTATE_CEN, Attribute(0));
- setAttribute(SAMPLE_OFF, Attribute(0));
- setAttribute(SAMPLE_X, Attribute(0.0));
- setAttribute(SAMPLE_Y, Attribute(0.0));
- setAttribute(SAMPLE_Z, Attribute(0.0));
- init();
+const int SCDPanelErrors::defaultIndexValue = 0;
+
+/// Constructor
+SCDPanelErrors::SCDPanelErrors() : m_setupFinished(false) {
+ declareParameter("XShift", 0.0, "Shift factor in X");
+ declareParameter("YShift", 0.0, "Shift factor in Y");
+ declareParameter("ZShift", 0.0, "Shift factor in Z");
+ declareParameter("XRotate", 0.0, "Rotate angle in X");
+ declareParameter("YRotate", 0.0, "Rotate angle in Y");
+ declareParameter("ZRotate", 0.0, "Rotate angle in Z");
+ declareAttribute("FileName", Attribute("", true));
+ declareAttribute("Workspace", Attribute(""));
+ declareAttribute("Bank", Attribute(""));
}
-
-void SCDPanelErrors::init() {
-
- declareParameter("f0_detWidthScale", 1.0, "panel Width");
- declareParameter("f0_detHeightScale", 1.0, "panel Height");
-
- declareParameter("f0_Xoffset", 0.0, "Panel Center x offset");
- declareParameter("f0_Yoffset", 0.0, "Panel Center y offset");
- declareParameter("f0_Zoffset", 0.0, "Panel Center z offset");
-
- declareParameter("f0_Xrot", 0.0,
- "Rotation(degrees) Panel Center in x axis direction");
- declareParameter("f0_Yrot", 0.0,
- "Rotation(degrees) Panel Center in y axis direction");
- declareParameter("f0_Zrot", 0.0,
- "Rotation(degrees) Panel Center in z axis direction");
-
- declareParameter("l0", 0.0, "Initial Flight Path");
- declareParameter("t0", 0.0, "Time offset");
- declareParameter("SampleX", 0.0, "Sample x offset");
- declareParameter("SampleY", 0.0, "Sample y offset");
- declareParameter("SampleZ", 0.0, "Sample z offset");
+/**
+ * The movedetector function changes detector position and angles
+ * @param x :: The shift along the X-axis
+ * @param y :: The shift along the Y-axis
+ * @param z :: The shift along the Z-axis
+ * @param rotx :: The rotation around the X-axis
+ * @param roty :: The rotation around the Y-axis
+ * @param rotz :: The rotation around the Z-axis
+ * @param detname :: The detector name
+ * @param inputW :: The workspace
+ */
+
+void SCDPanelErrors::moveDetector(
+ double x, double y, double z, double rotx, double roty, double rotz,
+ std::string detname, Workspace_sptr inputW) const {
+
+ IAlgorithm_sptr alg1 = Mantid::API::AlgorithmFactory::Instance().create("MoveInstrumentComponent", -1);
+ alg1->initialize();
+ alg1->setChild(true);
+ alg1->setLogging(false);
+ alg1->setProperty<Workspace_sptr>("Workspace", inputW);
+ alg1->setPropertyValue("ComponentName", detname);
+ // Move in cm for small shifts
+ alg1->setProperty("X", x * 0.01);
+ alg1->setProperty("Y", y * 0.01);
+ alg1->setProperty("Z", z * 0.01);
+ alg1->setPropertyValue("RelativePosition", "1");
+ alg1->execute();
+
+ IAlgorithm_sptr algx = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
+ algx->initialize();
+ algx->setChild(true);
+ algx->setLogging(false);
+ algx->setProperty<Workspace_sptr>("Workspace", inputW);
+ algx->setPropertyValue("ComponentName", detname);
+ algx->setProperty("X", 1.0);
+ algx->setProperty("Y", 0.0);
+ algx->setProperty("Z", 0.0);
+ algx->setProperty("Angle", rotx);
+ algx->setPropertyValue("RelativeRotation", "1");
+ algx->execute();
+
+ IAlgorithm_sptr algy = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
+ algy->initialize();
+ algy->setChild(true);
+ algy->setLogging(false);
+ algy->setProperty<Workspace_sptr>("Workspace", inputW);
+ algy->setPropertyValue("ComponentName", detname);
+ algy->setProperty("X", 0.0);
+ algy->setProperty("Y", 1.0);
+ algy->setProperty("Z", 0.0);
+ algy->setProperty("Angle", roty);
+ algy->setPropertyValue("RelativeRotation", "1");
+ algy->execute();
+
+ IAlgorithm_sptr algz = Mantid::API::AlgorithmFactory::Instance().create("RotateInstrumentComponent", -1);
+ algz->initialize();
+ algz->setChild(true);
+ algz->setLogging(false);
+ algz->setProperty<Workspace_sptr>("Workspace", inputW);
+ algz->setPropertyValue("ComponentName", detname);
+ algz->setProperty("X", 0.0);
+ algz->setProperty("Y", 0.0);
+ algz->setProperty("Z", 1.0);
+ algz->setProperty("Angle", rotz);
+ algz->setPropertyValue("RelativeRotation", "1");
+ algz->execute();
}
-void SCDPanelErrors::getPeaks() const {
- // if the pointer is set just return
- if (m_peaks && m_peaks->rowCount() > 0)
+/// Evaluate the function for a list of arguments and given scaling factor
+void SCDPanelErrors::eval(double xrotate, double yrotate, double zrotate, double xshift, double yshift, double zshift,
+ double *out, const double *xValues,
+ const size_t nData) const {
+ UNUSED_ARG(xValues);
+ if (nData == 0)
return;
- if (PeakName.empty())
- throw std::invalid_argument(
- "Cannot retrieve peaks workspace from empty string");
-
- // get the workspace from the framework
- m_peaks =
- AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(PeakName);
-
- // error check it
- if (!m_peaks || m_peaks->rowCount() < 1)
- throw std::invalid_argument("There are no peaks in the peaks workspace or "
- "no PeakWorkspace named \"" +
- PeakName + "\"");
-}
-
-void SCDPanelErrors::Check(DataObjects::PeaksWorkspace_sptr &pkwsp,
- const double *xValues, const size_t nData,
- size_t &StartX, size_t &EndX) const {
- // TODO need to error check this
- // if (NLatticeParametersSet < (int) nAttributes()-2)
- // {
- // g_log.error("Not all lattice parameters have been set");
- // throw std::invalid_argument("Not all lattice parameters have been set");
- // }
-
- if (!pkwsp) {
- throw std::invalid_argument("Cannot find a PeaksWorkspace ");
- }
-
- if (pkwsp->getNumberPeaks() < 4) {
- throw std::invalid_argument("Not enough peaks to fit ");
- }
-
- if ((m_startX > static_cast<int>(nData) - 1) ||
- (m_endX > static_cast<int>(nData) - 1)) {
- throw std::invalid_argument(X_START + " and " + X_END +
- " attributes are out of range");
- }
-
- StartX = 0;
- if (m_startX > 0)
- StartX = static_cast<size_t>(m_startX);
- EndX = nData - 1;
- if (m_endX > static_cast<int>(StartX))
- EndX = static_cast<size_t>(m_endX);
-
- if (xValues[StartX] != floor(xValues[StartX])) {
- throw std::invalid_argument("Improper workspace. xVals must be integer");
- }
-
- if (xValues[StartX] < 0 || xValues[StartX] >= pkwsp->rowCount()) {
-
- throw std::invalid_argument("Improper workspace. xVals correspond to an "
- "index in the PeaksWorkspace");
- }
-
- if ((EndX - StartX + 1) / 3 < 4) {
- throw std::invalid_argument("Not enough peaks to process banks " +
- BankNames);
+ setupData();
+
+ moveDetector(xrotate, yrotate, zrotate, xshift, yshift, zshift,m_bank,m_workspace);
+ DataObjects::PeaksWorkspace_sptr inputP =
+ boost::dynamic_pointer_cast<DataObjects::PeaksWorkspace>(m_workspace);
+ Geometry::Instrument_sptr inst = boost::const_pointer_cast<Geometry::Instrument>(inputP->getInstrument());
+ int j = 0;
+ for (size_t i = 0; i < nData; i += 3) {
+ DataObjects::Peak peak = inputP->getPeak(j);
+ j++;
+ V3D Q1 = peak.getQSampleFrame();
+ DataObjects::Peak peak2(inst,Q1,peak.getGoniometerMatrix());
+ peak2.setDetectorID(peak.getDetectorID());
+ V3D Q3 = peak2.getQSampleFrame();
+ out[i] = Q3[0];
+ out [i+1] = Q3[1];
+ out[i+2] = Q3[2];
}
+ moveDetector(-xrotate, -yrotate, -zrotate, -xshift, -yshift, -zshift,m_bank,m_workspace);
}
-Instrument_sptr
-SCDPanelErrors::getNewInstrument(const Geometry::IPeak &peak) const {
-
- Geometry::Instrument_const_sptr instSave = peak.getInstrument();
- auto pmap = boost::make_shared<ParameterMap>();
- boost::shared_ptr<const Geometry::ParameterMap> pmapSv =
- instSave->getParameterMap();
-
- if (!instSave) {
- g_log.error(" Not all peaks have an instrument");
- throw std::invalid_argument(" Not all peaks have an instrument");
- }
- auto instChange = boost::make_shared<Geometry::Instrument>();
-
- if (!instSave->isParametrized()) {
- boost::shared_ptr<Geometry::Instrument> instClone(instSave->clone());
- auto Pinsta = boost::make_shared<Geometry::Instrument>(instSave, pmap);
-
- instChange = Pinsta;
- } else // catch(...)
- {
- // TODO eliminate next 2 lines. not used
- boost::shared_ptr<const IComponent> inst3 =
- boost::dynamic_pointer_cast<const IComponent>(instSave);
- // updateParams(pmapSv, pmap, inst3);
-
- auto P1 = boost::make_shared<Geometry::Instrument>(
- instSave->baseInstrument(), pmap);
- instChange = P1;
- }
-
- if (!instChange) {
- g_log.error("Cannot 'clone' instrument");
- throw logic_error("Cannot clone instrument");
- }
- std::vector<std::string> GroupBanks;
-
- boost::split(GroupBanks, BankNames, boost::is_any_of("!"));
-
- for (size_t group = 0; group < (size_t)GroupBanks.size(); ++group) {
- string prefix = "f" + std::to_string(group) + "_";
-
- std::vector<std::string> bankNames;
- Quat rot = Quat(getParameter(prefix + "Xrot"), Kernel::V3D(1.0, 0.0, 0.0)) *
- Quat(getParameter(prefix + "Yrot"), Kernel::V3D(0.0, 1.0, 0.0)) *
- Quat(getParameter(prefix + "Zrot"), Kernel::V3D(0.0, 0.0, 1.0));
-
- boost::split(bankNames, GroupBanks[group], boost::is_any_of("/"));
-
- SCDCalibratePanels::FixUpBankParameterMap(
- bankNames, instChange,
- V3D(getParameter(prefix + "Xoffset"), getParameter(prefix + "Yoffset"),
- getParameter(prefix + "Zoffset")),
- rot, getParameter(prefix + "detWidthScale"),
- getParameter(prefix + "detHeightScale"), pmapSv, RotateCenters);
-
- } // for each group
-
- V3D SampPos = instChange->getSample()->getPos();
- SampPos[0] += getParameter("SampleX") + SampleX;
- SampPos[1] += getParameter("SampleY") + SampleY;
- SampPos[2] += getParameter("SampleZ") + SampleZ;
-
- SCDCalibratePanels::FixUpSourceParameterMap(instChange, getParameter("l0"),
- SampPos, pmapSv);
-
- return instChange;
+/**
+ * Calculate the function values.
+ * @param out :: The output buffer for the calculated values.
+ * @param xValues :: The array of x-values.
+ * @param nData :: The size of the data.
+ */
+void SCDPanelErrors::function1D(double *out, const double *xValues,
+ const size_t nData) const {
+ const double xrotate = getParameter("XRotate");
+ const double yrotate = getParameter("YRotate");
+ const double zrotate = getParameter("ZRotate");
+ const double xshift = getParameter("XShift");
+ const double yshift = getParameter("YShift");
+ const double zshift = getParameter("ZShift");
+ eval(xrotate, yrotate, zrotate, xshift, yshift, zshift, out, xValues, nData);
}
-Peak SCDPanelErrors::createNewPeak(const Geometry::IPeak &peak_old,
- Geometry::Instrument_sptr instrNew,
- double T0, double L0) {
- Geometry::Instrument_const_sptr inst = peak_old.getInstrument();
- if (inst->getComponentID() != instrNew->getComponentID()) {
- g_log.error("All peaks must have the same instrument");
- throw invalid_argument("All peaks must have the same instrument");
- }
-
- double T = peak_old.getTOF() + T0;
-
- int ID = peak_old.getDetectorID();
-
- Kernel::V3D hkl = peak_old.getHKL();
- // peak_old.setDetectorID(ID); //set det positions
- Peak peak(instrNew, ID, peak_old.getWavelength(), hkl,
- peak_old.getGoniometerMatrix());
-
- Wavelength wl;
-
- wl.initialize(L0, peak.getL2(), peak.getScattering(), 0,
- peak_old.getInitialEnergy(), 0.0);
-
- peak.setWavelength(wl.singleFromTOF(T));
- peak.setIntensity(peak_old.getIntensity());
- peak.setSigmaIntensity(peak_old.getSigmaIntensity());
- peak.setRunNumber(peak_old.getRunNumber());
- peak.setBinCount(peak_old.getBinCount());
-
- //!!!peak.setDetectorID(ID);
- return peak;
+/**
+ * function derivatives
+ * @param out :: The output Jacobian matrix: function derivatives over its
+ * parameters.
+ * @param xValues :: The function arguments
+ * @param nData :: The size of xValues.
+ */
+void SCDPanelErrors::functionDeriv1D(API::Jacobian *out,
+ const double *xValues,
+ const size_t nData) {
+ FunctionDomain1DView domain(xValues, nData);
+ this->calNumericalDeriv(domain, *out);
}
-void SCDPanelErrors::function1D(double *out, const double *xValues,
- const size_t nData) const {
- g_log.debug() << "Start function 1D\n";
-
- // return early if there is nothing to do
- if (nData == 0)
- return;
-
- if (!m_unitCell)
- throw runtime_error(
- "Cannot evaluate function without setting the lattice constants");
-
- // error check the parameters
- double r = checkForNonsenseParameters();
- if (r != 0) {
- for (size_t i = 0; i < nData; ++i)
- out[i] = 100 + r;
-
- g_log.debug() << "Parametersxx for " << BankNames << ">=";
- for (size_t i = 0; i < nParams(); ++i)
- g_log.debug() << getParameter(i) << ",";
- g_log.debug() << "\n";
-
- return;
- }
-
- // determine the range of data to fit by index
- size_t StartX;
- size_t EndX;
- this->getPeaks();
- Check(m_peaks, xValues, nData, StartX, EndX);
-
- g_log.debug() << "BankNames " << BankNames << " Number of peaks"
- << (EndX - StartX + 1) / 3 << '\n';
-
- // some pointers for the updated instrument
- boost::shared_ptr<Geometry::Instrument> instChange =
- getNewInstrument(m_peaks->getPeak(0));
-
- //---------------------------- Calculate q and hkl vectors-----------------
-
- vector<Kernel::V3D> hkl_vectors;
- vector<Kernel::V3D> q_vectors;
- double t0 = getParameter("t0");
- double l0 = getParameter("l0");
- for (size_t i = StartX; i <= EndX; i += 3) {
- // the x-values are the peak indices as triplets, convert them to size_t
- if (xValues[i] < 0.)
- throw invalid_argument(
- "Improper workspace. xVals must be positive integers");
- size_t pkIndex = std::lround(xValues[i]); // just round to nearest int
- if (pkIndex >= m_peaks->rowCount()) {
+/// Clear all data
+void SCDPanelErrors::clear() const {
+ m_setupFinished = false;
+}
- g_log.error() << "Improper workspace set " << pkIndex << "\n";
- throw invalid_argument("Improper workspace. xVals correspond to an index "
- "in the PeaksWorkspace");
+/** Set a value to attribute attName
+ * @param attName :: The attribute name
+ * @param value :: The new value
+ */
+void SCDPanelErrors::setAttribute(const std::string &attName,
+ const IFunction::Attribute &value) {
+ if (attName == "FileName") {
+ std::string fileName = value.asUnquotedString();
+ if (fileName.empty()) {
+ storeAttributeValue("FileName", Attribute("", true));
+ return;
}
-
- IPeak &peak_old = m_peaks->getPeak(static_cast<int>(pkIndex));
- Kernel::V3D hkl = peak_old.getHKL();
-
- // eliminate tolerance cause only those peaks that are OK should be here
- if (IndexingUtils::ValidIndex(hkl, tolerance)) {
- hkl_vectors.push_back(hkl);
- Peak peak = createNewPeak(peak_old, instChange, t0, l0);
- q_vectors.push_back(peak.getQSampleFrame());
+ FileValidator fval;
+ std::string error = fval.isValid(fileName);
+ if (error == "") {
+ storeAttributeValue(attName, Attribute(fileName, true));
+ storeAttributeValue("Workspace", Attribute(""));
+ } else {
+ // file not found
+ throw Kernel::Exception::FileError(error, fileName);
+ }
+ load(fileName);
+ } else if (attName == "Workspace") {
+ std::string wsName = value.asString();
+ if (!wsName.empty()) {
+ storeAttributeValue(attName, value);
+ storeAttributeValue("FileName", Attribute("", true));
+ loadWorkspace(wsName);
}
+ } else {
+ IFunction::setAttribute(attName, value);
+ m_setupFinished = false;
}
+}
- //----------------------------------Calculate out
- //----------------------------------
-
- // determine the OrientedLattice for converting to Q-sample
- Geometry::OrientedLattice lattice(*m_unitCell);
- lattice.setUB(m_peaks->sample().getOrientedLattice().getUB());
-
- // cumulative error
- double chiSq = 0; // for debug log message
-
- for (size_t i = 0; i < StartX; ++i)
- out[i] = 0.;
- for (size_t i = 0; i < q_vectors.size(); ++i) {
- Kernel::V3D err = q_vectors[i] - lattice.qFromHKL(hkl_vectors[i]);
- size_t outIndex = 3 * i + StartX;
- out[outIndex + 0] = err[0];
- out[outIndex + 1] = err[1];
- out[outIndex + 2] = err[2];
- chiSq += err[0] * err[0] + err[1] * err[1] + err[2] * err[2];
+/**
+ * Load input file as a Nexus file.
+ * @param fname :: The file name
+ */
+void SCDPanelErrors::load(const std::string &fname) {
+ IAlgorithm_sptr loadAlg =
+ Mantid::API::AlgorithmFactory::Instance().create("Load", -1);
+ loadAlg->initialize();
+ loadAlg->setChild(true);
+ loadAlg->setLogging(false);
+ try {
+ loadAlg->setPropertyValue("Filename", fname);
+ loadAlg->setPropertyValue("OutputWorkspace",
+ "_SCDPanelErrors_fit_data_");
+ loadAlg->execute();
+ } catch (std::runtime_error &) {
+ throw std::runtime_error(
+ "Unable to load Nexus file for SCDPanelErrors function.");
}
- for (size_t i = EndX; i < nData; ++i)
- out[i] = 0.;
-
- g_log.debug() << "Parameters\n";
-
- for (size_t i = 0; i < this->nParams(); ++i)
- g_log.debug() << setw(20) << parameterName(i) << setw(20) << getParameter(i)
- << '\n';
-
- g_log.debug() << " chi Squared=" << std::setprecision(12) << chiSq
- << '\n';
-
- // Get values for test program. TODO eliminate
- g_log.debug() << " out[evenxx]=";
- for (size_t i = 0; i < std::min<size_t>(nData, 30); ++i)
- g_log.debug() << out[i] << " ";
-
- g_log.debug() << '\n';
+ Workspace_sptr ws = loadAlg->getProperty("OutputWorkspace");
+ API::Workspace_sptr resData =
+ boost::dynamic_pointer_cast<Mantid::API::Workspace>(ws);
+ loadWorkspace(resData);
}
-double SCDPanelErrors::checkForNonsenseParameters() const {
-
- double Dwdth = getParameter(0);
- double Dhght = getParameter(1);
- double x = getParameter(2);
- double y = getParameter(3);
- double z = getParameter(4);
- double rx = getParameter(5);
- double ry = getParameter(6);
- double rz = getParameter(7);
- double L0 = getParameter(8);
- double T0 = getParameter(9);
-
- double r = 0.;
- if (L0 < 1.)
- r = 1. - L0;
-
- if (fabs(T0) > 20.)
- r += (T0 - 20.) * 2.;
-
- if (Dwdth < .5 || Dwdth > 2.)
- r += 3. * fabs(1 - Dwdth);
-
- if (Dhght < .5 || Dhght > 2.)
- r += 3. * fabs(1. - Dhght);
-
- if (fabs(x) > .35)
- r += fabs(x) * .2;
-
- if (fabs(y) > .35)
- r += fabs(y) * .2;
-
- if (fabs(z) > .35)
- r += fabs(z) * .2;
-
- if (fabs(rx) > 15.)
- r += fabs(rx) * .02;
-
- if (fabs(ry) > 15.)
- r += fabs(ry) * .02;
-
- if (fabs(rz) > 15.)
- r += fabs(rz) * .02;
-
- return 5. * r;
+/**
+ * Load the points from a PeaksWorkspace
+ * @param wsName :: The workspace to load from
+ */
+void SCDPanelErrors::loadWorkspace(const std::string &wsName) const {
+ auto ws = AnalysisDataService::Instance().retrieveWS<API::Workspace>(wsName);
+ loadWorkspace(ws);
}
-void updateDerivResult(PeaksWorkspace_sptr peaks, V3D &unRotDeriv,
- Matrix<double> &Result, size_t peak,
- vector<int> &peakIndx) {
- Matrix<double> GonMatrix =
- peaks->getPeak(peakIndx[peak]).getGoniometerMatrix();
- GonMatrix.Invert();
-
- V3D RotDeriv = GonMatrix * unRotDeriv;
-
- for (int kk = 0; kk < 3; ++kk)
- Result[kk][peak] = RotDeriv[kk];
+/**
+ * Load the points from a PeaksWorkspace
+ * @param ws :: The workspace to load from
+ */
+void SCDPanelErrors::loadWorkspace(
+ boost::shared_ptr<API::Workspace> ws) const {
+ m_workspace = ws;
+ m_setupFinished = false;
}
-void SCDPanelErrors::functionDeriv1D(Jacobian *out, const double *xValues,
- const size_t nData) {
- if (nData <= 0)
+/**
+ * Fill in the workspace name and bank
+ */
+void SCDPanelErrors::setupData() const {
+ if (m_setupFinished) {
+ g_log.debug() << "Re-setting isn't required.";
return;
- FunctionDomain1DView domain(xValues, nData);
- calNumericalDeriv(domain, *out);
-}
-
-DataObjects::Workspace2D_sptr
-SCDPanelErrors::calcWorkspace(DataObjects::PeaksWorkspace_sptr &pwks,
- std::vector<std::string> &bankNames,
- double tolerance) {
- int N = 0;
- if (tolerance < 0)
- tolerance = .5;
- tolerance = std::min<double>(.5, tolerance);
-
- Mantid::MantidVec xRef;
-
- for (auto &bankName : bankNames)
- for (size_t j = 0; j < pwks->rowCount(); ++j) {
- Geometry::IPeak &peak = pwks->getPeak(static_cast<int>(j));
- if (peak.getBankName().compare(bankName) == 0)
- if (peak.getH() != 0 || peak.getK() != 0 || peak.getL() != 0)
- if (peak.getH() - floor(peak.getH()) < tolerance ||
- floor(peak.getH() + 1) - peak.getH() < tolerance)
- if (peak.getK() - floor(peak.getK()) < tolerance ||
- floor(peak.getK() + 1) - peak.getK() < tolerance)
- if (peak.getL() - floor(peak.getL()) < tolerance ||
- floor(peak.getL() + 1) - peak.getL() < tolerance) {
- N++;
- xRef.push_back(static_cast<double>(j));
- xRef.push_back(static_cast<double>(j));
- xRef.push_back(static_cast<double>(j));
- }
- }
-
- MatrixWorkspace_sptr mwkspc = API::WorkspaceFactory::Instance().create(
- "Workspace2D", static_cast<size_t>(3), 3 * N, 3 * N);
-
- auto pX = Kernel::make_cow<HistogramData::HistogramX>(std::move(xRef));
- mwkspc->setX(0, pX);
- mwkspc->setX(1, pX);
- mwkspc->setX(2, pX);
-
- return boost::dynamic_pointer_cast<DataObjects::Workspace2D>(mwkspc);
-}
-
-void SCDPanelErrors::setAttribute(const std::string &attName,
- const Attribute &value) {
- if (!hasAttribute(attName))
- throw std::invalid_argument("Not a valid attribute name \"" + attName +
- "\"");
-
- bool recalcB = false;
- if (attName == LATTICE_A) {
- a = value.asDouble();
- a_set = true;
- recalcB = true;
- } else if (attName == LATTICE_B) {
- b = value.asDouble();
- b_set = true;
- recalcB = true;
- } else if (attName == LATTICE_C) {
- c = value.asDouble();
- c_set = true;
- recalcB = true;
- } else if (attName == LATTICE_ALPHA) {
- alpha = value.asDouble();
- alpha_set = true;
- recalcB = true;
- } else if (attName == LATTICE_BETA) {
- beta = value.asDouble();
- beta_set = true;
- recalcB = true;
- } else if (attName == LATTICE_GAMMA) {
- gamma = value.asDouble();
- gamma_set = true;
- recalcB = true;
- } else if (attName == (BANK_NAMES)) {
- BankNames = value.asString();
- BankNames_set = true;
- } else if (attName == PEAKS_WKSP) {
- PeakName = value.asString();
- PeakName_set = true;
- } else if (attName == NUM_GROUPS) {
- if (NGroups_set) {
- g_log.error("Cannot set NGroups more than once");
- throw new std::invalid_argument("Cannot set NGroups more than once");
- }
- NGroups = value.asInt();
- for (int k = 1; k < NGroups; ++k) {
- std::string prefix = "f" + std::to_string(k) + "_";
- declareParameter(prefix + "detWidthScale", 1.0, "panel Width");
- declareParameter(prefix + "detHeightScale", 1.0, "panelHeight");
-
- declareParameter(prefix + "Xoffset", 0.0, "Panel Center x offset");
- declareParameter(prefix + "Yoffset", 0.0, "Panel Center y offset");
- declareParameter(prefix + "Zoffset", 0.0, "Panel Center z offset");
-
- declareParameter(prefix + "Xrot", 0.0,
- "Rotation(degrees) Panel Center in x axis direction");
- declareParameter(prefix + "Yrot", 0.0,
- "Rotation(degrees) Panel Center in y axis direction");
- declareParameter(prefix + "Zrot", 0.0,
- "Rotation(degrees) Panel Center in z axis direction");
- }
- NGroups_set = true;
- } else if (attName == ROTATE_CEN) {
- int v = value.asInt();
- if (v == 0)
- RotateCenters = false;
- else
- RotateCenters = true;
+ }
- } else if (attName == SAMPLE_OFF) {
- int v = value.asInt();
- if (v == 0)
- SampleOffsets = false;
+ if (!m_workspace) {
+ std::string wsName = getAttribute("Workspace").asString();
+ if (wsName.empty())
+ throw std::invalid_argument("Data not set for function " + this->name());
else
- SampleOffsets = true;
+ loadWorkspace(wsName);
+ }
- } else if (attName == SAMPLE_X) {
- SampleX = value.asDouble();
- sampleX_set = true;
- } else if (attName == SAMPLE_Y) {
- SampleY = value.asDouble();
- sampleY_set = true;
- } else if (attName == SAMPLE_Z) {
- SampleZ = value.asDouble();
- sampleZ_set = true;
- } else if (attName == X_START) {
- m_startX = value.asInt();
- startX_set = true;
- } else if (attName == "endX") {
- m_endX = value.asInt();
- endX_set = true;
- } else
- throw std::invalid_argument("Not a valid attribute name \"" + attName +
- "\"");
+ m_bank = getAttribute("Bank").asString();
- if (recalcB) {
- if (a_set && b_set && c_set && alpha_set && beta_set && gamma_set) {
- m_unitCell =
- boost::make_shared<Geometry::UnitCell>(a, b, c, alpha, beta, gamma);
- }
- }
-}
+ g_log.debug() << "Setting up " << m_workspace->name() << " bank " << m_bank
+ << '\n';
-boost::shared_ptr<const Geometry::IComponent>
-SCDPanelErrors::findBank(std::string bankName) {
- return bankDetMap.find(bankName)->second;
+ m_setupFinished = true;
}
} // namespace Crystal
diff --git a/Framework/DataHandling/src/MoveInstrumentComponent.cpp b/Framework/DataHandling/src/MoveInstrumentComponent.cpp
index ef67846d94c..8f7f7ecaf64 100644
--- a/Framework/DataHandling/src/MoveInstrumentComponent.cpp
+++ b/Framework/DataHandling/src/MoveInstrumentComponent.cpp
@@ -5,6 +5,7 @@
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/Exception.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"
+#include "MantidDataObjects/PeaksWorkspace.h"
namespace Mantid {
namespace DataHandling {
@@ -23,7 +24,7 @@ MoveInstrumentComponent::MoveInstrumentComponent() {}
void MoveInstrumentComponent::init() {
// When used as a Child Algorithm the workspace name is not used - hence the
// "Anonymous" to satisfy the validator
- declareProperty(make_unique<WorkspaceProperty<MatrixWorkspace>>(
+ declareProperty(make_unique<WorkspaceProperty<Workspace>>(
"Workspace", "Anonymous", Direction::InOut),
"The name of the workspace for which the new instrument "
"configuration will have an effect. Any other workspaces "
@@ -50,8 +51,32 @@ void MoveInstrumentComponent::init() {
* @throw std::runtime_error Thrown with Workspace problems
*/
void MoveInstrumentComponent::exec() {
- // Get the workspace
- MatrixWorkspace_sptr WS = getProperty("Workspace");
+ // Get the input workspace
+ Workspace_sptr ws = getProperty("Workspace");
+ MatrixWorkspace_sptr inputW =
+ boost::dynamic_pointer_cast<MatrixWorkspace>(ws);
+ DataObjects::PeaksWorkspace_sptr inputP =
+ boost::dynamic_pointer_cast<DataObjects::PeaksWorkspace>(ws);
+
+ // Get some stuff from the input workspace
+ Instrument_sptr inst;
+ if (inputW) {
+ inst = boost::const_pointer_cast<Instrument>(inputW->getInstrument());
+ if (!inst)
+ throw std::runtime_error("Could not get a valid instrument from the "
+ "MatrixWorkspace provided as input");
+ } else if (inputP) {
+ inst = boost::const_pointer_cast<Instrument>(inputP->getInstrument());
+ if (!inst)
+ throw std::runtime_error("Could not get a valid instrument from the "
+ "PeaksWorkspace provided as input");
+ } else {
+ if (!inst)
+ throw std::runtime_error("Could not get a valid instrument from the "
+ "workspace and it does not seem to be valid as "
+ "input (must be either MatrixWorkspace or "
+ "PeaksWorkspace");
+ }
const std::string ComponentName = getProperty("ComponentName");
const int DetID = getProperty("DetectorID");
const double X = getProperty("X");
@@ -59,10 +84,6 @@ void MoveInstrumentComponent::exec() {
const double Z = getProperty("Z");
const bool relativePosition = getProperty("RelativePosition");
- // Get the ParameterMap reference before the instrument so that
- // we avoid a copy
- Geometry::ParameterMap &pmap = WS->instrumentParameters();
- Instrument_const_sptr inst = WS->getInstrument();
IComponent_const_sptr comp;
// Find the component to move
@@ -92,8 +113,15 @@ void MoveInstrumentComponent::exec() {
TransformType positionType = Absolute;
if (relativePosition)
positionType = Relative;
- Geometry::ComponentHelper::moveComponent(*comp, pmap, V3D(X, Y, Z),
- positionType);
+ if (inputW) {
+ Geometry::ParameterMap &pmap = inputW->instrumentParameters();
+ Geometry::ComponentHelper::moveComponent(*comp, pmap, V3D(X, Y, Z),
+ positionType);
+ } else if (inputP) {
+ Geometry::ParameterMap &pmap = inputP->instrumentParameters();
+ Geometry::ComponentHelper::moveComponent(*comp, pmap, V3D(X, Y, Z),
+ positionType);
+ }
}
} // namespace DataHandling
--
GitLab