Commit ffd01776 authored by Peterson, Peter's avatar Peterson, Peter
Browse files

Fix behavior for negative time-of-flight to d

Apparently some instruments still have a DAS/DAQ/DAE that sometimes
creates a negative time-of-flight in the raw data file. This
demonstrated that the quadratic equation was not correctly converting to
d-spacing when DIFA==0. This PR adds an extra trap to simplify the
calculation in this case. Also changes the order of special cases to
move calculating intermediate results until as late as possible. In
principle, this should execute the case of difa==0 faster as zero is not
squared then square rooted.
parent 718b953a
......@@ -644,56 +644,80 @@ void dSpacing::init() {
}
}
}
// force the assumption that difc is positive in debug builds
assert(difc > 0.);
}
double dSpacing::singleToTOF(const double x) const {
if (!isInitialized())
throw std::runtime_error("dSpacingBase::singleToTOF called before object "
"has been initialized.");
return difa * x * x + difc * x + tzero;
if (difa == 0.)
return difc * x + tzero;
else
return difa * x * x + difc * x + tzero;
}
/**
* DIFA * d^2 + DIFC * d + T0 - TOF = 0
*
* Use the citardauq formula to solve quadratic in order to minimise loss of precision. citarqauq (quadratic spelled
* backwards) is an alternate formulation of the quadratic formula. DIFC and sqrt term are often similar and the
* "classic" quadratic formula involves calculating their difference in the numerator
*
* 2*(T0 - TOF) (T0 - TOF)
* d = ------------------------------------------- = ---------------------------------------------------
* -DIFC -+ SQRT(DIFC^2 - 4*DIFA*(T0 - TOF)) 0.5 * DIFC (-1 -+ SQRT(1 - 4*DIFA*(T0 - TOF)/DIFC^2)
*
* the variables in this formulation are the same as the quadratic formula
* a = difa square term
* b = DIFC linear term - assumed to be positive
* c = T0 - TOF constant term
*/
double dSpacing::singleFromTOF(const double tof) const {
// DIFA * d^2 + DIFC * d + T0 - TOF = 0
// Use the citardauq formula to solve quadratic in order to minimise loss of precision. DIFC and sqrt term are often
// similar and the "classic" quadratic formula involves calculating their difference in the numerator
// 2*(T0 - TOF) (T0 - TOF)
// d = ------------------------------------------- = ---------------------------------------------------
// -DIFC -+ SQRT(DIFC^2 - 4*DIFA*(T0 - TOF)) 0.5 * DIFC (-1 -+ SQRT(1 - 4*DIFA*(T0 - TOF)/DIFC^2)
// dealing with various edge cases
if (!isInitialized())
throw std::runtime_error("dSpacingBase::singleFromTOF called before object "
"has been initialized.");
if (!toDSpacingError.empty())
throw std::runtime_error(toDSpacingError);
double c = tzero - tof;
// force the assumption that difc is positive in debug builds
assert(difc > 0.);
// handle special cases first...
// citardauq formula hides non-zero root if c=0
// don't need to solve a quadratic when difa==0
if (difa == 0.)
return (tof - tzero) / difc;
// citardauq formula hides non-zero root if tof==tzero
if (tof == tzero) {
if (difa < 0)
if (difa < 0.)
return -difc / difa;
else
return 0;
return 0.;
}
double sqrtTerm = 1 - 4 * difa * c / (difc * difc);
if (sqrtTerm < 0) {
throw std::runtime_error("Cannot convert to d spacing. Quadratic doesn't have real roots");
}
if ((difa > 0) && (c > 0)) {
// non-physical result
if ((difa > 0.) && (tzero > tof)) {
throw std::runtime_error("Cannot convert to d spacing. Quadratic doesn't "
"have a positive root");
}
// general equation
const double linearTerm = tzero - tof;
const double sqrtTerm = 1 - 4 * difa * linearTerm / (difc * difc);
if (sqrtTerm < 0.) {
throw std::runtime_error("Cannot convert to d spacing. Quadratic doesn't have real roots");
}
// pick smallest positive root. Since difc is positive it just depends on sign of c
// Note - c is generally negative
if (c > 0)
if (linearTerm > 0)
// single positive root
return c / (0.5 * difc * (-1 + sqrt(sqrtTerm)));
return linearTerm / (0.5 * difc * (-1 + sqrt(sqrtTerm)));
else
// two positive roots. pick most negative denominator to get smallest root
return c / (0.5 * difc * (-1 - sqrt(sqrtTerm)));
return linearTerm / (0.5 * difc * (-1 - sqrt(sqrtTerm)));
}
double dSpacing::conversionTOFMin() const {
......
......@@ -577,13 +577,19 @@ public:
}
void testdSpacing_fromTOF() {
std::vector<double> x(1, 1001.1), y(1, 1.0);
std::vector<double> yy = y;
const std::vector<double> x_in{-1., 0., 1001.1, 16000.};
const std::vector<double> y_in{1., 2., 3., 4.};
std::vector<double> x(x_in.begin(), x_in.end());
std::vector<double> y(y_in.begin(), y_in.end());
double difc =
2.0 * Mantid::PhysicalConstants::NeutronMass * sin(0.5) * (1.0 + 1.0) * 1e-4 / Mantid::PhysicalConstants::h;
TS_ASSERT_THROWS_NOTHING(d.fromTOF(x, y, 1.0, 1, {{UnitParams::difc, difc}}))
TS_ASSERT_DELTA(x[0], 2.065172, 0.000001)
TS_ASSERT(yy == y)
TS_ASSERT_THROWS_NOTHING(d.fromTOF(x, y, 1.0, 1, {{UnitParams::difc, difc}}));
TS_ASSERT(y == y_in);
for (size_t i = 0; i < x.size(); ++i)
TS_ASSERT_DELTA(x[i], x_in[i] / difc, 0.000001);
// test for exception thrown
x[0] = 1.0;
TS_ASSERT_THROWS(d.fromTOF(x, y, 1.0, 1, {{UnitParams::difc, -1.0}}), const std::runtime_error &)
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment