Commit 1c81861a authored by Gigg, Martyn Anthony's avatar Gigg, Martyn Anthony
Browse files

Improve efficiency in handling goniometer. Refs #5397

Added reference returns for goniometer & its axes plus moved where
the computation of the matrix occurs.
parent ee812b46
......@@ -116,16 +116,15 @@ namespace Mantid
/// Returns the bin boundaries for a given value
std::pair<double, double> histogramBinBoundaries(const double energyValue) const;
/** @return a reference to the Goniometer object for this run */
Mantid::Geometry::Goniometer & getGoniometer()
{ return m_goniometer; }
/** @return a reference to the const Goniometer object for this run */
const Mantid::Geometry::Goniometer & getGoniometer() const
{ return m_goniometer; }
/// Set the gonoimeter & read the values from the logs if told to do so
void setGoniometer(const Geometry::Goniometer & goniometer, const bool useLogValues);
/** @return A reference to the const Goniometer object for this run */
inline const Geometry::Goniometer & getGoniometer() const { return m_goniometer; }
/** @return A reference to the non-const Goniometer object for this run */
inline Geometry::Goniometer & mutableGoniometer() { return m_goniometer; }
// Retrieve the goniometer rotation matrix
Mantid::Kernel::DblMatrix getGoniometerMatrix();
const Kernel::DblMatrix & getGoniometerMatrix() const;
/**
* Add a log entry
......@@ -164,6 +163,8 @@ namespace Mantid
private:
/// Adds all the time series in from one property manager into another
void mergeMergables(Mantid::Kernel::PropertyManager & sum, const Mantid::Kernel::PropertyManager & toAdd);
/// Calculate the gonoimeter matrix
void calculateGoniometerMatrix();
/// Static reference to the logger class
static Kernel::Logger &g_log;
......
......@@ -168,33 +168,6 @@ Kernel::Logger& Run::g_log = Kernel::Logger::get("Run");
}
}
/** Adds all the time series in the second property manager to those in the first
* @param sum the properties to add to
* @param toAdd the properties to add
*/
void Run::mergeMergables(Mantid::Kernel::PropertyManager & sum, const Mantid::Kernel::PropertyManager & toAdd)
{
// get pointers to all the properties on the right-handside and prepare to loop through them
const std::vector<Property*> inc = toAdd.getProperties();
std::vector<Property*>::const_iterator end = inc.end();
for (std::vector<Property*>::const_iterator it=inc.begin(); it != end;++it)
{
const std::string rhs_name = (*it)->name();
try
{
//now get pointers to the same properties on the left-handside
Property * lhs_prop(sum.getProperty(rhs_name));
lhs_prop->merge(*it);
}
catch (Exception::NotFoundError &)
{
//copy any properties that aren't already on the left hand side
Property * copy = (*it)->clone();
//And we add a copy of that property to *this
sum.declareProperty(copy, "");
}
}
}
//-----------------------------------------------------------------------------------------------
/**
......@@ -556,6 +529,26 @@ Kernel::Logger& Run::g_log = Kernel::Logger::get("Run");
}
//-----------------------------------------------------------------------------------------------
/**
* Set the gonoimeter & optionally read the values from the logs
* @param goniometer :: A refernce to a goniometer
* @param useLogValues :: If true, recalculate the goniometer using the log values
*/
void Run::setGoniometer(const Geometry::Goniometer & goniometer, const bool useLogValues)
{
Geometry::Goniometer old = m_goniometer;
try
{
m_goniometer = goniometer; //copy it in
if(useLogValues) calculateGoniometerMatrix();
}
catch(std::runtime_error&)
{
m_goniometer = old;
throw;
}
}
/** Get the gonimeter rotation matrix, calculated using the
* previously set Goniometer object as well as the angles
* loaded in the run (if any).
......@@ -564,49 +557,11 @@ Kernel::Logger& Run::g_log = Kernel::Logger::get("Run");
*
* @return 3x3 double rotation matrix
*/
Mantid::Kernel::DblMatrix Run::getGoniometerMatrix()
const Mantid::Kernel::DblMatrix & Run::getGoniometerMatrix() const
{
for (size_t i=0; i < m_goniometer.getNumberAxes(); ++i)
{
std::string name = m_goniometer.getAxis(i).name;
if (this->hasProperty(name))
{
Property * prop = this->getProperty(name);
//try time series property
TimeSeriesProperty<double> * tsp = dynamic_cast<TimeSeriesProperty<double> *>(prop);
if (tsp)
{
// Set that angle
m_goniometer.setRotationAngle(i, tsp->getStatistics().mean);
g_log.debug() << "Goniometer angle " << name << " set to " <<tsp->getStatistics().mean << std::endl;
}
else
{
//try Property with value
PropertyWithValue<double> *pvp=dynamic_cast<PropertyWithValue<double> *>(prop);
if (pvp)
{
// Set that angle
m_goniometer.setRotationAngle(i, boost::lexical_cast<double>(pvp->value()));
g_log.debug() << "Goniometer angle " << name << " set to " << boost::lexical_cast<double>(pvp->value()) << std::endl;
}
else
{
throw std::runtime_error("Sample log for goniometer angle '" + name + "' was not a TimeSeriesProperty<double> or PropertyWithValue<double>.");
}
}
}
else
throw std::runtime_error("Could not find goniometer angle '" + name + "' in the run sample logs.");
}
return m_goniometer.getR();
}
//--------------------------------------------------------------------------------------------
/** Save the object to an open NeXus file.
* @param file :: open NeXus file
......@@ -706,6 +661,52 @@ Kernel::Logger& Run::g_log = Kernel::Logger::get("Run");
}
}
//-----------------------------------------------------------------------------------------------------------------------
// Private methods
//-----------------------------------------------------------------------------------------------------------------------
/** Adds all the time series in the second property manager to those in the first
* @param sum the properties to add to
* @param toAdd the properties to add
*/
void Run::mergeMergables(Mantid::Kernel::PropertyManager & sum, const Mantid::Kernel::PropertyManager & toAdd)
{
// get pointers to all the properties on the right-handside and prepare to loop through them
const std::vector<Property*> inc = toAdd.getProperties();
std::vector<Property*>::const_iterator end = inc.end();
for (std::vector<Property*>::const_iterator it=inc.begin(); it != end;++it)
{
const std::string rhs_name = (*it)->name();
try
{
//now get pointers to the same properties on the left-handside
Property * lhs_prop(sum.getProperty(rhs_name));
lhs_prop->merge(*it);
}
catch (Exception::NotFoundError &)
{
//copy any properties that aren't already on the left hand side
Property * copy = (*it)->clone();
//And we add a copy of that property to *this
sum.declareProperty(copy, "");
}
}
}
/**
* Calculate the goniometer matrix
*/
void Run::calculateGoniometerMatrix()
{
for (size_t i=0; i < m_goniometer.getNumberAxes(); ++i)
{
const double angle = getLogAsSingleValue(m_goniometer.getAxis(i).name, Kernel::Math::Mean);
m_goniometer.setRotationAngle(i, angle);
}
}
/// @cond
/// Macro to instantiate concrete template members
#define INSTANTIATE(TYPE) \
......
......@@ -377,7 +377,9 @@ public:
// No axes by default
TS_ASSERT_EQUALS( runInfo.getGoniometer().getNumberAxes(), 0 );
// Now does copy work?
runInfo.getGoniometer().makeUniversalGoniometer();
Goniometer gm;
gm.makeUniversalGoniometer();
runInfo.setGoniometer(gm, false);
TS_ASSERT_EQUALS( runInfo.getGoniometer().getNumberAxes(), 3 );
Run runCopy(runInfo);
TS_ASSERT_EQUALS( runCopy.getGoniometer().getNumberAxes(), 3 );
......@@ -394,6 +396,30 @@ public:
runInfo.addProperty(tsp);
}
/** Setting up a goniometer and the angles to feed it
* using sample logs, then getting the right rotation matrix out.
*/
void test_setGoniometerWhenLogsDoNotExistsThrows()
{
Run runInfo;
Goniometer gm;
gm.makeUniversalGoniometer();
TS_ASSERT_THROWS(runInfo.setGoniometer(gm, true), std::runtime_error);
}
/** Setting up a goniometer and the angles to feed it
* using sample logs, then getting the right rotation matrix out.
*/
void test_setGoniometer_Not_Using_Logs_Preserves_Input()
{
Run runInfo;
DblMatrix rotation(3,3, true);
Goniometer gm(rotation);
TS_ASSERT_EQUALS(runInfo.getGoniometer().getNumberAxes(), 0);
TS_ASSERT_EQUALS(runInfo.getGoniometer().getR(), rotation);
}
/** Setting up a goniometer and the angles to feed it
* using sample logs, then getting the right rotation matrix out.
......@@ -404,7 +430,9 @@ public:
addTimeSeriesEntry(runInfo, "phi", 90.0);
addTimeSeriesEntry(runInfo, "chi", 90.0);
addTimeSeriesEntry(runInfo, "omega", 90.0);
runInfo.getGoniometer().makeUniversalGoniometer();
Goniometer gm;
gm.makeUniversalGoniometer();
runInfo.setGoniometer(gm, true);
DblMatrix r = runInfo.getGoniometerMatrix();
V3D rot = r * V3D(-1,0,0);
TS_ASSERT_EQUALS( rot, V3D(1, 0, 0));
......@@ -418,7 +446,10 @@ public:
addTimeSeriesEntry(runInfo, "phi", 45.0);
addTimeSeriesEntry(runInfo, "chi", 90.0);
addTimeSeriesEntry(runInfo, "omega", 0.0);
runInfo.getGoniometer().makeUniversalGoniometer();
Goniometer gm;
gm.makeUniversalGoniometer();
runInfo.setGoniometer(gm, true);
DblMatrix r = runInfo.getGoniometerMatrix();
V3D rot = r * V3D(-1,0,0);
TS_ASSERT_EQUALS( rot, V3D(0, -sqrt(0.5), sqrt(0.5)));
......@@ -432,7 +463,6 @@ public:
th.createFile("RunTest.nxs");
Run run1;
run1.getGoniometer().makeUniversalGoniometer();
addTimeSeriesEntry(run1, "double_series", 45.0);
run1.addProperty( new PropertyWithValue<int>("int_val", 1234) );
run1.addProperty( new PropertyWithValue<std::string>("string_val", "help_im_stuck_in_a_log_file") );
......@@ -442,6 +472,10 @@ public:
addTimeSeriesEntry(run1, "omega", 78.9);
addTimeSeriesEntry(run1, "proton_charge", 78.9);
Goniometer gm;
gm.makeUniversalGoniometer();
run1.setGoniometer(gm, true);
run1.storeHistogramBinBoundaries(m_test_energy_bins);
run1.saveNexus(th.file, "logs");
......
......@@ -125,10 +125,8 @@ namespace Crystal
if (gon.getNumberAxes() == 0)
g_log.warning() << "Empty goniometer created; will always return an identity rotation matrix." << std::endl;
// All went well, copy the goniometer into it
ws->mutableRun().getGoniometer() = gon;
//force it to read the values
ws->mutableRun().getGoniometerMatrix(); //it will throw if log values are not found
// All went well, copy the goniometer into it. It will throw if the log values cannot be found
ws->mutableRun().setGoniometer(gon, true);
}
......
......@@ -67,7 +67,7 @@ public:
TS_ASSERT( alg.isExecuted() ); //no log values
// Check the results
Goniometer & gon = ws->mutableRun().getGoniometer();
const Goniometer & gon = ws->mutableRun().getGoniometer();
TS_ASSERT_EQUALS( gon.getNumberAxes(), 0);
DblMatrix rot = ws->mutableRun().getGoniometerMatrix();
TSM_ASSERT_EQUALS( "Goniometer Rotation matrix is 3x3 identity", rot, DblMatrix(3,3, true) );
......@@ -95,7 +95,7 @@ public:
TS_ASSERT( alg.isExecuted() ); //no log values
// Check the results
Goniometer & gon = ws->mutableRun().getGoniometer();
const Goniometer & gon = ws->mutableRun().getGoniometer();
TS_ASSERT_EQUALS( gon.getNumberAxes(), 2);
TS_ASSERT_EQUALS( gon.getAxis(0).name, "angle1");
......@@ -134,7 +134,7 @@ public:
else
{
// Check the results
Goniometer & gon = ws->mutableRun().getGoniometer();
const Goniometer & gon = ws->mutableRun().getGoniometer();
TS_ASSERT_EQUALS( gon.getNumberAxes(), numExpected);
}
AnalysisDataService::Instance().remove("SetGoniometerTest_ws");
......
......@@ -1457,8 +1457,6 @@ EventWorkspace_sptr LoadEventNexus::createEmptyEventWorkspace()
eventWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
eventWS->setYUnit("Counts");
// Create a default "Universal" goniometer in the Run object
eventWS->mutableRun().getGoniometer().makeUniversalGoniometer();
return eventWS;
}
......@@ -1676,7 +1674,19 @@ BankPulseTimes * LoadEventNexus::runLoadNexusLogs(const std::string &nexusfilena
localWorkspace->mutableRun().addProperty("run_start", run_start.toISO8601String(), true );
}
else
{
alg->getLogger().warning() << "Empty proton_charge sample log. You will not be able to filter by time.\n";
}
/// Attempt to make a gonoimeter from the logs
try
{
Goniometer gm;
gm.makeUniversalGoniometer();
localWorkspace->mutableRun().setGoniometer(gm, true);
}
catch(std::runtime_error &)
{
}
}
catch (...)
{
......
......@@ -308,10 +308,11 @@ namespace DataHandling
outputWS->mutableRun().addLogData(new PropertyWithValue<std::string>("ki_over_kf_scaling", kikfscaling==1?"true":"false"));
//Set Goniometer
Geometry::Goniometer gm;
gm.pushAxis("psi",0,1,0,psi);
outputWS->mutableRun().setGoniometer(gm, true);
outputWS->mutableRun().getGoniometer().pushAxis("psi",0,1,0,psi);
//generate instrument
Geometry::Instrument_sptr instrument(new Geometry::Instrument("NXSPE"));
outputWS->setInstrument(instrument);
......
......@@ -110,10 +110,9 @@ namespace Geometry
// Set rotation angle for an axis
void setRotationAngle( size_t axisnumber, double value);
// Get axis object
GoniometerAxis getAxis(size_t axisnumber);
const GoniometerAxis getAxis(size_t axisnumber) const;
const GoniometerAxis & getAxis(size_t axisnumber) const;
// Get axis object
GoniometerAxis getAxis(std::string axisname);
const GoniometerAxis & getAxis(std::string axisname) const;
// Return the number of axes
size_t getNumberAxes() const;
// Make a default universal goniometer
......
......@@ -159,26 +159,17 @@ void Goniometer::setRotationAngle( size_t axisnumber, double value)
/// Get GoniometerAxis obfject using motor number
/// @param axisnumber :: axis number (from 0)
GoniometerAxis Goniometer::getAxis( size_t axisnumber)
{
if (axisnumber >= motors.size())
throw std::out_of_range("Goniometer::getAxis(): axis number specified is too large.");
return motors.at(axisnumber);//it will throw out of range exception if axisnumber is not in range
}
/// Get GoniometerAxis obfject using motor number
/// @param axisnumber :: axis number (from 0)
const GoniometerAxis Goniometer::getAxis( size_t axisnumber) const
const GoniometerAxis & Goniometer::getAxis( size_t axisnumber) const
{
return motors.at(axisnumber);//it will throw out of range exception if axisnumber is not in range
}
/// Get GoniometerAxis object using motor name
/// @param axisname :: axis name
GoniometerAxis Goniometer::getAxis( std::string axisname)
const GoniometerAxis & Goniometer::getAxis( std::string axisname) const
{
bool found=false;
std::vector<GoniometerAxis>::iterator it;
for(it=motors.begin(); it<motors.end(); ++it)
for(auto it=motors.begin(); it<motors.end(); ++it)
{
if(axisname.compare((*it).name)==0)
{
......
......@@ -21,9 +21,8 @@
File change history is stored at: <https://svn.mantidproject.org/mantid/trunk/Code/Mantid>.
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
#include "MantidAPI/ParamFunctionAttributeHolder.h"
#include "MantidAPI/ExperimentInfo.h"
#include "MantidAPI/ParamFunctionAttributeHolder.h"
#include "MantidMDAlgorithms/Quantification/ForegroundModelFactory.h"
namespace Mantid
......@@ -35,9 +34,9 @@ namespace Mantid
* to be combined with a resolution calculation.
*
* A concrete model should override the following functions
* - declareParamters() : Defines the parameters within the model
* - declareParameters() : Defines the parameters within the model
* - scatteringIntensity() : Returns a value for the cross section with the
* current parameters
* current parameters
*/
class DLLExport ForegroundModel: public API::ParamFunctionAttributeHolder
{
......@@ -50,9 +49,13 @@ namespace Mantid
/// Constructor taking the fitted function to access the current parameter values
ForegroundModel(const API::IFunction & fittingFunction);
/// Returns the type of model
virtual ModelType modelType() const = 0;
/// Calculates the intensity for the model for the current parameters, expt description & ND point
virtual double scatteringIntensity(const API::ExperimentInfo & exptSetup, const std::vector<double> & point) const = 0;
/// Set a reference to the convolved fitting function. Needed as we need a default constructor
void setFunctionUnderMinimization(const API::IFunction & fitFunction);
/// Declares the parameters
void declareParameters();
/// Return the initial value of the parameter according to the fit by index
......@@ -64,12 +67,6 @@ namespace Mantid
/// Return the current parameter according to the fit by name
double getCurrentParameterValue(const std::string& name) const;
/// Returns the type of model
virtual ModelType modelType() const = 0;
/// Calculates the intensity for the model for the current parameters, expt description & ND point
virtual double scatteringIntensity(const API::ExperimentInfo & exptSetup, const std::vector<double> & point) const = 0;
protected:
/// Returns a reference to the fitting function
const API::IFunction & functionUnderMinimization() const;
......@@ -81,12 +78,10 @@ namespace Mantid
void function(const Mantid::API::FunctionDomain&, Mantid::API::FunctionValues&) const {}
/// Hide these
using ParamFunctionAttributeHolder::getParameter;
using ParamFunction::getParameter;
/// A (non-owning) pointer to the function undergoing fitting
const API::IFunction * m_fittingFunction;
/// A (non-owning) pointer to an experiment info
const API::ExperimentInfo * m_exptInfo;
/// An offset for the number of parameters that were declared before this one
size_t m_parOffset;
};
......
......@@ -45,14 +45,21 @@ namespace Mantid
std::string name() const { return "Strontium122"; }
/// Declare the fitting parameters
void declareParameters();
/// Declare fixed attributes
void declareAttributes();
/// Called when an attribute is set
void setAttribute(const std::string & name, const API::IFunction::Attribute& attr);
/// Returns the type of model
ModelType modelType() const { return Broad; }
/// Calculates the intensity for the model for the current parameters.
double scatteringIntensity(const API::ExperimentInfo & exptDescr, const std::vector<double> & point) const;
/// Twin type attribute
int m_twinType;
/// MultEps attribute
bool m_multEps;
/// Magnetic form factor cache
PhysicalConstants::MagneticFormFactorTable m_formFactorTable;
};
}
......
......@@ -21,7 +21,7 @@ namespace Mantid
* @param A reference to the fitting function
*/
ForegroundModel::ForegroundModel(const API::IFunction & fittingFunction)
: m_fittingFunction(NULL), m_parOffset(0)
: API::ParamFunctionAttributeHolder(), m_fittingFunction(NULL), m_parOffset(0)
{
setFunctionUnderMinimization(fittingFunction);
}
......
......@@ -20,11 +20,16 @@ namespace Mantid
namespace // anonymous
{
/// Enumerate parameter positions
enum { Seff = 0, J1a = 1, J1b = 2, J2 = 3, SJc = 4, GammaSlope = 5, MultEps = 6, TwinType = 7};
enum { Seff = 0, J1a = 1, J1b = 2, J2 = 3, SJc = 4, GammaSlope = 5 };
/// Number of parameters
const unsigned int NPARAMS = 8;
const unsigned int NPARAMS = 6;
/// Parameter names, same order as above
const char * PAR_NAMES[8] = { "Seff", "J1a", "J1b", "J2", "SJc", "GammaSlope", "MultEps", "TwinType" };
const char * PAR_NAMES[NPARAMS] = { "Seff", "J1a", "J1b", "J2", "SJc", "GammaSlope" };
/// N attrs
const unsigned int NATTS = 2;
/// Attribute names
const char * ATTR_NAMES[NATTS] = { "MultEps", "TwinType" };
}
Strontium122::Strontium122()
......@@ -42,6 +47,35 @@ namespace Mantid
}
}
/**
* Declare fixed attributes
*/
void Strontium122::declareAttributes()
{
for(unsigned int i = 0; i < NATTS; ++i)
{
declareAttribute(ATTR_NAMES[i], API::IFunction::Attribute(1));
}
}
/**
* Called when an attribute is set
* @param name :: The name of the attribute
* @param attr :: The value of the attribute
*/
void Strontium122::setAttribute(const std::string & name, const API::IFunction::Attribute& attr)
{
int asInt = attr.asInt();
if(name == ATTR_NAMES[0])
{
m_multEps = (asInt > 0);
}
else if(name == ATTR_NAMES[1])
{
m_twinType = asInt;
}
}
/**
* Calculates the scattering intensity
* @param exptSetup :: Details of the current experiment
......@@ -54,13 +88,29 @@ namespace Mantid
const double qsqr = qx*qx + qy*qy + qz*qz;
const double epssqr = eps*eps;
const Geometry::OrientedLattice & latticeRot = exptSetup.sample().getOrientedLattice();
static const double twoPi = 2.*M_PI;
const Geometry::OrientedLattice & lattice = exptSetup.sample().getOrientedLattice();
const Kernel::DblMatrix & gr = exptSetup.run().getGoniometerMatrix();
const Kernel::DblMatrix & ub = lattice.getUB();
const Kernel::DblMatrix & ub = latticeRot.getUB();
const double ub00(ub[0][0]), ub01(ub[0][1]), ub02(ub[0][2]),
ub10(ub[1][0]), ub11(ub[1][1]), ub12(ub[1][2]),
ub20(ub[2][0]), ub21(ub[2][1]), ub22(ub[2][2]);
double ub00(0.0), ub01(0.0), ub02(0.0),
ub10(0.0), ub11(0.0), ub12(0.0),
ub20(0.0), ub21(0.0), ub22(0.0);
for(unsigned int i = 0; i < 3; ++i)
{
ub00 += gr[0][i]*ub[i][0];
ub01 += gr[0][i]*ub[i][1];
ub02 += gr[0][i]*ub[i][2];
ub10 += gr[1][i]*ub[i][0];
ub11 += gr[