Commit 2191c259 authored by Danny Hindson's avatar Danny Hindson Committed by Zhang, Chen
Browse files

Extra way of initialising dSpacing and MomentumTransfer

Enable initialising dSpacing and MomentumTransfer using L2 and two theta
as an alternative to supplying the diff constants. This might be more
user friendly and it also moves the difc formula into the Unit class

Add utility functions to manage extra params map
parent 78937a43
......@@ -70,11 +70,10 @@ public:
double azimuthal(const size_t index) const;
std::pair<double, double> geographicalAngles(const size_t index) const;
Kernel::V3D position(const size_t index) const;
std::tuple<double, double, double>
Kernel::UnitParametersMap
diffractometerConstants(const size_t index,
std::vector<detid_t> &uncalibratedDets) const;
std::tuple<double, double, double>
diffractometerConstants(const size_t index) const;
Kernel::UnitParametersMap diffractometerConstants(const size_t index) const;
double difcUncalibrated(const size_t index) const;
bool hasDetectors(const size_t index) const;
bool hasUniqueDetector(const size_t index) const;
......
......@@ -145,9 +145,9 @@ Kernel::V3D SpectrumInfo::position(const size_t index) const {
* @param warningDets A vector containing the det ids where an uncalibrated
* value was used in the situation where some dets have calibrated values and
* some don't
* @return tuple containing the average constants
* @return map containing the average constants
*/
std::tuple<double, double, double>
UnitParametersMap
SpectrumInfo::diffractometerConstants(const size_t index,
std::vector<detid_t> &warningDets) const {
if (m_detectorInfo.isScanning()) {
......@@ -177,37 +177,40 @@ SpectrumInfo::diffractometerConstants(const size_t index,
// if no calibration is found then return difc only based on the average
// of the detector L2 and twoThetas.
if (calibratedDets.size() == 0) {
return {0., difcUncalibrated(index), 0.};
return {{UnitParams::difc, difcUncalibrated(index)}};
}
return {difa / static_cast<double>(spectrumDefinition(index).size()),
difc / static_cast<double>(spectrumDefinition(index).size()),
tzero / static_cast<double>(spectrumDefinition(index).size())};
return {{UnitParams::difa,
difa / static_cast<double>(spectrumDefinition(index).size())},
{UnitParams::difc,
difc / static_cast<double>(spectrumDefinition(index).size())},
{UnitParams::tzero,
tzero / static_cast<double>(spectrumDefinition(index).size())}};
}
/** Calculate average diffractometer constants (DIFA, DIFC, TZERO) of detectors
* associated with this spectrum. Use calibrated values where possible, filling
* in with uncalibrated values where they're missing
/** Calculate average diffractometer constants (DIFA, DIFC, TZERO) of
* detectors associated with this spectrum. Use calibrated values where
* possible, filling in with uncalibrated values where they're missing
* @param index Index of the spectrum that constants are required for
* @return tuple containing the average constants
* @return map containing the average constants
*/
std::tuple<double, double, double>
UnitParametersMap
SpectrumInfo::diffractometerConstants(const size_t index) const {
std::vector<int> warningDets;
return diffractometerConstants(index, warningDets);
}
/** Calculate average uncalibrated DIFC value of detectors associated with this
* spectrum
/** Calculate average uncalibrated DIFC value of detectors associated with
* this spectrum
* @param index Index of the spectrum that DIFC is required for
* @return The average DIFC
*/
double SpectrumInfo::difcUncalibrated(const size_t index) const {
// calculate difc based on the average of the detector L2 and twoThetas. This
// will be different to the average of the per detector difcs. This is for
// backwards compatibility because Mantid always used to calculate spectrum
// level difc's this way
return 1. / Mantid::Geometry::Conversion::tofToDSpacingFactor(
l1(), l2(index), twoTheta(index), 0.);
// calculate difc based on the average of the detector L2 and twoThetas.
// This will be different to the average of the per detector difcs. This is
// for backwards compatibility because Mantid always used to calculate
// spectrum level difc's this way
return 1. / Kernel::Units::tofToDSpacingFactor(l1(), l2(index),
twoTheta(index), 0.);
}
/** Get the detector values relevant to unit conversion for a workspace index
......@@ -217,9 +220,9 @@ double SpectrumInfo::difcUncalibrated(const size_t index) const {
* @param signedTheta :: Return twotheta with sign or without
* @param wsIndex :: The workspace index
* @param pmap :: a map containing values for conversion parameters that are
required by unit classes to perform their conversions eg efixed. It can contain
values on the way in if a look up isn't desired here eg if value supplied in
parameters to the calling algorithm
required by unit classes to perform their conversions eg efixed. It can
contain values on the way in if a look up isn't desired here eg if value
supplied in parameters to the calling algorithm
*/
void SpectrumInfo::getDetectorValues(const Kernel::Unit &inputUnit,
const Kernel::Unit &outputUnit,
......@@ -239,7 +242,6 @@ void SpectrumInfo::getDetectorValues(const Kernel::Unit &inputUnit,
pmap[UnitParams::twoTheta] = twoTheta(wsIndex);
} catch (const std::runtime_error &e) {
g_log.warning(e.what());
pmap[UnitParams::twoTheta] = std::numeric_limits<double>::quiet_NaN();
}
if (emode != Kernel::DeltaEMode::Elastic &&
pmap.find(UnitParams::efixed) == pmap.end()) {
......@@ -257,17 +259,13 @@ void SpectrumInfo::getDetectorValues(const Kernel::Unit &inputUnit,
std::vector<detid_t> warnDetIds;
try {
std::vector<std::string> diffConstUnits = {"dSpacing", "MomentumTransfer",
"Empty"};
std::set<std::string> diffConstUnits = {"dSpacing", "MomentumTransfer",
"Empty"};
if ((emode == Kernel::DeltaEMode::Elastic) &&
(std::find(diffConstUnits.begin(), diffConstUnits.end(),
inputUnit.unitID()) != diffConstUnits.end()) &&
(std::find(diffConstUnits.begin(), diffConstUnits.end(),
outputUnit.unitID()) != diffConstUnits.end())) {
auto [difa, difc, tzero] = diffractometerConstants(wsIndex, warnDetIds);
pmap[UnitParams::difa] = difa;
pmap[UnitParams::difc] = difc;
pmap[UnitParams::tzero] = tzero;
(diffConstUnits.count(inputUnit.unitID()) ||
diffConstUnits.count(outputUnit.unitID()))) {
auto diffConstsMap = diffractometerConstants(wsIndex, warnDetIds);
pmap.insert(diffConstsMap.begin(), diffConstsMap.end());
if (warnDetIds.size() > 0) {
createDetectorIdLogMessages(warnDetIds, wsIndex);
}
......@@ -307,7 +305,8 @@ void SpectrumInfo::createDetectorIdLogMessages(
detIDstring);
}
/// Returns true if the spectrum is associated with detectors in the instrument.
/// Returns true if the spectrum is associated with detectors in the
/// instrument.
bool SpectrumInfo::hasDetectors(const size_t index) const {
// Workspaces can contain invalid detector IDs. Those IDs will be silently
// ignored here until this is fixed.
......
......@@ -300,11 +300,12 @@ void populateTable(ITableWorkspace_sptr &t, const MatrixWorkspace_sptr &ws,
if (isMonitor) {
colValues << 0. << 0. << 0. << 0.;
} else {
auto [difaValue, difcValue, tzeroValue] =
spectrumInfo.diffractometerConstants(wsIndex);
auto diffConsts = spectrumInfo.diffractometerConstants(wsIndex);
auto difcValueUncalibrated = spectrumInfo.difcUncalibrated(wsIndex);
colValues << difaValue << difcValue << difcValueUncalibrated
<< tzeroValue;
// map will create an entry with zero value if not present already
colValues << diffConsts[UnitParams::difa]
<< diffConsts[UnitParams::difc] << difcValueUncalibrated
<< diffConsts[UnitParams::tzero];
}
}
} catch (const std::exception &) {
......
......@@ -328,11 +328,10 @@ void SaveGSS::generateBankHeader(std::stringstream &out, const API::SpectrumInfo
const auto l1 = spectrumInfo.l1();
const auto l2 = spectrumInfo.l2(specIndex);
const auto twoTheta = spectrumInfo.twoTheta(specIndex);
double difc;
std::tie(std::ignore, difc, std::ignore) =
spectrumInfo.diffractometerConstants(specIndex);
auto diffConstants = spectrumInfo.diffractometerConstants(specIndex);
out << "# Total flight path " << (l1 + l2) << "m, tth "
<< (twoTheta * 180. / M_PI) << "deg, DIFC " << difc << "\n";
<< (twoTheta * 180. / M_PI) << "deg, DIFC "
<< diffConstants[Kernel::UnitParams::difc] << "\n";
}
out << "# Data for spectrum :" << specIndex << "\n";
}
......
......@@ -20,6 +20,7 @@
#include "MantidKernel/Exception.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Unit.h"
#include <memory>
#include <nexus/NeXusFile.hpp>
......@@ -1290,20 +1291,9 @@ namespace Conversion {
* @param offset
* @return
*/
double tofToDSpacingFactor(const double l1, const double l2, const double twoTheta, const double offset) {
if (offset <= -1.) // not physically possible, means result is negative d-spacing
{
std::stringstream msg;
msg << "Encountered offset of " << offset << " which converts data to negative d-spacing\n";
throw std::logic_error(msg.str());
}
auto sinTheta = std::sin(twoTheta / 2);
const double numerator = (1.0 + offset);
sinTheta *= (l1 + l2);
return (numerator * CONSTANT) / sinTheta;
double tofToDSpacingFactor(const double l1, const double l2,
const double twoTheta, const double offset) {
return Kernel::Units::tofToDSpacingFactor(l1, l2, twoTheta, offset);
}
} // namespace Conversion
......
......@@ -15,6 +15,7 @@
#include "MantidKernel/EigenConversionHelpers.h"
#include "MantidKernel/Exception.h"
#include "MantidKernel/MultiThreaded.h"
#include "MantidKernel/Unit.h"
namespace Mantid {
namespace Geometry {
......@@ -295,8 +296,8 @@ std::tuple<double, double, double> DetectorInfo::diffractometerConstants(
}
double DetectorInfo::difcUncalibrated(const size_t index) const {
return 1. / Mantid::Geometry::Conversion::tofToDSpacingFactor(
l1(), l2(index), twoTheta(index), 0.);
return 1. / Kernel::Units::tofToDSpacingFactor(l1(), l2(index),
twoTheta(index), 0.);
}
std::pair<double, double>
......
......@@ -25,6 +25,8 @@ namespace Kernel {
enum class UnitParams { l2, twoTheta, efixed, delta, difa, difc, tzero };
// list of parameter names and values - some code may relay on map behaviour
// where [] creates element with value 0 if param name not present
using UnitParametersMap = std::unordered_map<UnitParams, double>;
/** The base units (abstract) class. All concrete units should inherit from
......@@ -232,12 +234,12 @@ protected:
bool initialized;
/// l1 :: The source-sample distance (in metres)
double l1;
/// l2 :: The sample-detector distance (in metres)
double l2;
/// emode :: The energy mode (0=elastic, 1=direct geometry, 2=indirect
/// geometry)
int emode;
/// additional parameters
/// l2 :: distance from sample to detector (in metres)
/// twoTheta :: scattering angle in radians
/// efixed :: Value of fixed energy: EI (emode=1) or EF (emode=2) (in meV)
/// difc :: diffractometer constant DIFC
const UnitParametersMap *m_params;
......@@ -407,18 +409,15 @@ protected:
double factorFrom; ///< Constant factor for from conversion
};
//=================================================================================================
/// d-Spacing in Angstrom
class MANTID_KERNEL_DLL dSpacing : public Unit {
public:
const std::string unitID() const override; ///< "dSpacing"
const std::string caption() const override { return "d-Spacing"; }
const UnitLabel label() const override;
MANTID_KERNEL_DLL double tofToDSpacingFactor(const double l1, const double l2,
const double twoTheta,
const double offset);
class MANTID_KERNEL_DLL dSpacingBase : public Unit {
public:
double singleToTOF(const double x) const override;
double singleFromTOF(const double tof) const override;
void init() override;
Unit *clone() const override;
double conversionTOFMin() const override;
double conversionTOFMax() const override;
double calcTofMin(const double difc, const double difa, const double tzero,
......@@ -427,7 +426,7 @@ public:
const double tofmax = 0.);
/// Constructor
dSpacing();
dSpacingBase();
protected:
void validateUnitParams(const int emode,
......@@ -439,6 +438,20 @@ protected:
double valueThatGivesSmallestTOF;
};
//=================================================================================================
/// d-Spacing in Angstrom
class MANTID_KERNEL_DLL dSpacing : public dSpacingBase {
public:
const std::string unitID() const override; ///< "dSpacing"
const std::string caption() const override { return "d-Spacing"; }
const UnitLabel label() const override;
Unit *clone() const override;
/// Constructor
dSpacing();
};
//=================================================================================================
/// d-SpacingPerpendicular in Angstrom
class MANTID_KERNEL_DLL dSpacingPerpendicular : public Unit {
......@@ -469,7 +482,7 @@ protected:
//=================================================================================================
/// Momentum Transfer in Angstrom^-1
class MANTID_KERNEL_DLL MomentumTransfer : public dSpacing {
class MANTID_KERNEL_DLL MomentumTransfer : public dSpacingBase {
public:
const std::string unitID() const override; ///< "MomentumTransfer"
const std::string caption() const override { return "q"; }
......
......@@ -8,19 +8,41 @@
// Includes
//----------------------------------------------------------------------
#include "MantidKernel/Unit.h"
#include "MantidKernel/Logger.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/UnitLabelTypes.h"
#include <cfloat>
#include <sstream>
namespace Mantid {
namespace Kernel {
namespace {
// static logger object
Logger g_log("Unit");
bool ParamPresent(const UnitParametersMap &params, UnitParams param) {
return params.find(param) != params.end();
}
bool ParamPresentAndSet(const UnitParametersMap *params, UnitParams param,
double &var) {
auto it = params->find(param);
if (it != params->end()) {
var = it->second;
return true;
} else {
return false;
}
}
} // namespace
/**
* Default constructor
* Gives the unit an empty UnitLabel
*/
Unit::Unit() : initialized(false), l1(0), l2(0), emode(0) {}
Unit::Unit() : initialized(false), l1(0), emode(0) {}
bool Unit::operator==(const Unit &u) const { return unitID() == u.unitID(); }
......@@ -309,14 +331,12 @@ const UnitLabel Wavelength::label() const { return Symbol::Angstrom; }
void Wavelength::validateUnitParams(const int emode,
const UnitParametersMap &params) {
auto it = params.find(UnitParams::l2);
if (it == params.end()) {
if (!ParamPresent(params, UnitParams::l2)) {
throw std::runtime_error("An l2 value must be supplied in the extra "
"parameters when initialising " +
this->unitID() + " for conversion via TOF");
}
it = params.find(UnitParams::efixed);
if ((emode != 0) && (it == params.end())) {
if ((emode != 0) && (!ParamPresent(params, UnitParams::efixed))) {
throw std::runtime_error("An efixed value must be supplied in the extra "
"parameters when initialising " +
this->unitID() + " for conversion via TOF");
......@@ -325,19 +345,14 @@ void Wavelength::validateUnitParams(const int emode,
void Wavelength::init() {
// ------------ Factors to convert TO TOF ---------------------
double l2 = 0.0;
double ltot = 0.0;
double TOFisinMicroseconds = 1e6;
double toAngstroms = 1e10;
sfpTo = 0.0;
auto it = m_params->find(UnitParams::efixed);
if (it != m_params->end()) {
efixed = it->second;
}
it = m_params->find(UnitParams::l2);
if (it != m_params->end()) {
l2 = it->second;
}
ParamPresentAndSet(m_params, UnitParams::efixed, efixed);
ParamPresentAndSet(m_params, UnitParams::l2, l2);
if (emode == 1) {
ltot = l2;
......@@ -439,10 +454,8 @@ Energy::Energy() : Unit(), factorTo(DBL_MIN), factorFrom(DBL_MIN) {
}
void Energy::init() {
auto it = m_params->find(UnitParams::l2);
if (it != m_params->end()) {
l2 = it->second;
}
double l2 = 0.0;
ParamPresentAndSet(m_params, UnitParams::l2, l2);
{
const double TOFinMicroseconds = 1e6;
factorTo = sqrt(PhysicalConstants::NeutronMass / (2.0 * PhysicalConstants::meV)) * (l1 + l2) * TOFinMicroseconds;
......@@ -498,8 +511,7 @@ Energy_inWavenumber::Energy_inWavenumber() : Unit(), factorTo(DBL_MIN), factorFr
void Energy_inWavenumber::validateUnitParams(const int,
const UnitParametersMap &params) {
auto it = params.find(UnitParams::l2);
if (it == params.end()) {
if (!ParamPresent(params, UnitParams::l2)) {
throw std::runtime_error("An l2 value must be supplied in the extra "
"parameters when initialising " +
this->unitID() + " for conversion via TOF");
......@@ -507,10 +519,8 @@ void Energy_inWavenumber::validateUnitParams(const int,
}
void Energy_inWavenumber::init() {
auto it = m_params->find(UnitParams::l2);
if (it != m_params->end()) {
l2 = it->second;
}
double l2 = 0.0;
ParamPresentAndSet(m_params, UnitParams::l2, l2);
{
const double TOFinMicroseconds = 1e6;
factorTo =
......@@ -551,27 +561,51 @@ Unit *Energy_inWavenumber::clone() const { return new Energy_inWavenumber(*this)
*
* Conversion uses Bragg's Law: 2d sin(theta) = n * lambda
*/
DECLARE_UNIT(dSpacing)
const UnitLabel dSpacing::label() const { return Symbol::Angstrom; }
const double CONSTANT = (PhysicalConstants::h * 1e10) /
(2.0 * PhysicalConstants::NeutronMass * 1e6);
dSpacing::dSpacing()
: Unit(), difa(0), difc(DBL_MIN), tzero(0),
valueThatGivesLargestTOF(DBL_MAX), valueThatGivesSmallestTOF(0.) {
const double factor = 2.0 * M_PI;
addConversion("MomentumTransfer", factor, -1.0);
addConversion("QSquared", (factor * factor), -2.0);
/**
* Calculate and return conversion factor from tof to d-spacing.
* @param l1
* @param l2
* @param twoTheta scattering angle
* @param offset
* @return
*/
double tofToDSpacingFactor(const double l1, const double l2,
const double twoTheta, const double offset) {
if (offset <=
-1.) // not physically possible, means result is negative d-spacing
{
std::stringstream msg;
msg << "Encountered offset of " << offset
<< " which converts data to negative d-spacing\n";
throw std::logic_error(msg.str());
}
auto sinTheta = std::sin(twoTheta / 2);
const double numerator = (1.0 + offset);
sinTheta *= (l1 + l2);
return (numerator * CONSTANT) / sinTheta;
}
void dSpacing::validateUnitParams(const int, const UnitParametersMap &params) {
auto it = params.find(UnitParams::difc);
if (it == params.end()) {
throw std::runtime_error(
"A difc value must be supplied in the extra parameters when "
"initialising " +
this->unitID() + " for conversion via TOF");
}
if (it->second < 0.) {
dSpacingBase::dSpacingBase()
: Unit(), difa(0), difc(DBL_MIN), tzero(0),
valueThatGivesLargestTOF(DBL_MAX), valueThatGivesSmallestTOF(0.) {}
void dSpacingBase::validateUnitParams(const int,
const UnitParametersMap &params) {
double difc;
if (!ParamPresentAndSet(&params, UnitParams::difc, difc)) {
if (!ParamPresent(params, UnitParams::twoTheta) ||
(!ParamPresent(params, UnitParams::l2)))
throw std::runtime_error("A difc value or L2\two theta must be supplied "
"in the extra parameters when initialising " +
this->unitID() + " for conversion via TOF");
} else if (difc < 0.) {
throw std::runtime_error(
"A positive difc value must be supplied in the extra parameters when "
"initialising " +
......@@ -579,35 +613,44 @@ void dSpacing::validateUnitParams(const int, const UnitParametersMap &params) {
}
}
void dSpacing::init() {
void dSpacingBase::init() {
// First the crux of the conversion
difa = 0.;
difc = 0.;
tzero = 0.;
auto it = m_params->find(UnitParams::difc);
if (it != m_params->end()) {
difc = it->second;
}
it = m_params->find(UnitParams::difa);
if (it != m_params->end()) {
difa = it->second;
}
it = m_params->find(UnitParams::tzero);
if (it != m_params->end()) {
tzero = it->second;
ParamPresentAndSet(m_params, UnitParams::difa, difa);
ParamPresentAndSet(m_params, UnitParams::tzero, tzero);
if (!ParamPresentAndSet(m_params, UnitParams::difc, difc)) {
// also support inputs as L2, two theta
double l2;
if (ParamPresentAndSet(m_params, UnitParams::l2, l2)) {
double twoTheta;
if (ParamPresentAndSet(m_params, UnitParams::twoTheta, twoTheta)) {
if (difa != 0.) {
g_log.warning("Supplied difa ignored");
difa = 0.;
}
difc = 1. / tofToDSpacingFactor(l1, l2, twoTheta, 0.);
if (tzero != 0.) {
g_log.warning("Supplied tzero ignored");
tzero = 0.;
}
}
}
}
}
double dSpacing::singleToTOF(const double x) const {
double dSpacingBase::singleToTOF(const double x) const {
if (!isInitialized())
throw std::runtime_error(
"dSpacing::singleToTOF called before object has been initialized.");
throw std::runtime_error("dSpacingBase::singleToTOF called before object "
"has been initialized.");
return difa * x * x + difc * x + tzero;
}
double dSpacing::singleFromTOF(const double tof) const {
double dSpacingBase::singleFromTOF(const double tof) const {
if (!isInitialized())
throw std::runtime_error(
"dSpacing::singleFromTOF called before object has been initialized.");
throw std::runtime_error("dSpacingBase::singleFromTOF called before object "
"has been initialized.");
// handle special cases first...
if (tof == tzero) {
if (difa != 0)
......@@ -622,8 +665,8 @@ double dSpacing::singleFromTOF(const double tof) const {
"Cannot convert to d spacing. Quadratic doesn't have real roots");
}
if ((difa > 0) && ((tzero - tof) > 0)) {
throw std::runtime_error(
"Cannot convert to d spacing. Quadratic doesn't have a positive root");
throw std::runtime_error("Cannot convert to d spacing. Quadratic doesn't "
"have a positive root");
}
// and then solve quadratic using Muller formula
double sqrtTerm;
......@@ -638,7 +681,7 @@ double dSpacing::singleFromTOF(const double tof) const {
else
return (tof - tzero) / (0.5 * (difc + sqrtTerm));
}
double dSpacing::conversionTOFMin() const {
double dSpacingBase::conversionTOFMin() const {
// quadratic only has a min if difa is positive
if (difa > 0) {
// min of the quadratic is at d=-difc/(2*difa)
......@@ -652,7 +695,7 @@ double dSpacing::conversionTOFMin() const {
return TOFmin;
}
}
double dSpacing::conversionTOFMax() const {
double dSpacingBase::conversionTOFMax() const {
// quadratic only has a max if difa is negative
if (difa < 0) {
return std::min(DBL_MAX, tzero - difc * difc / (4 * difa));
......@@ -666,8 +709,8 @@ double dSpacing::conversionTOFMax() const {
}
}
double dSpacing::calcTofMin(const double difc, const double difa,
const double tzero, const double tofmin) {
double dSpacingBase::calcTofMin(const double difc, const double difa,
const double tzero, const double tofmin) {
Kernel::UnitParametersMap params{{Kernel::UnitParams::difa, difa},
{Kernel::UnitParams::difc, difc},
{Kernel::UnitParams::tzero, tzero}};
......@@ -675,8 +718,8 @@ double dSpacing::calcTofMin(const double difc, const double difa,
return std::max(conversionTOFMin(), tofmin);
}
double dSpacing::calcTofMax(const double difc, const double difa,
const double tzero, const double tofmax) {
double dSpacingBase::calcTofMax(const double difc, const double difa,
const double tzero, const double tofmax) {