Skip to content

Commit

Permalink
Refs #11847. Added unit test and modified LoadSpiceXML behavior.
Browse files Browse the repository at this point in the history
  • Loading branch information
@wdzhou
wdzhou committed May 28, 2015
1 parent ef55256 commit a2e2261
Show file tree
Hide file tree
Showing 6 changed files with 304 additions and 33 deletions.
1 change: 1 addition & 0 deletions Code/Mantid/Framework/API/src/ExperimentInfo.cpp
View file
Expand Up @@ -141,6 +141,7 @@ const std::string ExperimentInfo::toString() const {
*/
void ExperimentInfo::setInstrument(const Instrument_const_sptr &instr) {
if (instr->isParametrized()) {
g_log.notice("[DB] ExperimentInfo.setInstrument as being isParametrized.");
sptr_instrument = instr->baseInstrument();
m_parmap = instr->getParameterMap();
} else {
Expand Down
29 changes: 18 additions & 11 deletions Code/Mantid/Framework/DataHandling/src/LoadSpiceXML2DDet.cpp
View file
Expand Up @@ -265,6 +265,24 @@ void LoadSpiceXML2DDet::exec() {
if (spicetablewsname.size() > 0) {
setupSampleLogFromSpiceTable(outws, spicetablews, ptnumber);
}
else
{
// Set up 2theta from _2theta
if (outws->run().hasProperty("_2theta"))
{
Kernel::DateAndTime anytime(1000);
double logvalue = atof(outws->run().getProperty("_2theta")->value().c_str());
g_log.information() << "Set 2theta from _2thea with value " << logvalue << "\n";
TimeSeriesProperty<double> *newlogproperty =
new TimeSeriesProperty<double>("2theta");
newlogproperty->addValue(anytime, logvalue);
outws->mutableRun().addProperty(newlogproperty);
}
else
{
g_log.warning("No 2theta is set up for loading instrument.");
}
}

if (loadinstrument) {
loadInstrument(outws, idffilename);
Expand Down Expand Up @@ -613,20 +631,9 @@ bool LoadSpiceXML2DDet::getHB3AWavelength(MatrixWorkspace_sptr dataws,

void LoadSpiceXML2DDet::setXtoLabQ(API::MatrixWorkspace_sptr dataws,
const double &wavelength) {
// Geometry information
// Kernel::V3D sourcepos = dataws->getInstrument()->getSource()->getPos();
// Kernel::V3D samplepos = dataws->getInstrument()->getSample()->getPos();
// Kernel::V3D sample_source_vector = sourcepos - samplepos;

size_t numspec = dataws->getNumberHistograms();
for (size_t iws = 0; iws < numspec; ++iws) {
/*
Kernel::V3D detpos = dataws->getDetector(iws)->getPos();
Kernel::V3D detpos_sample_vector = detpos - samplepos;
double scattering_rad = detpos_sample_vector.angle(sample_source_vector);
*/

// double ki = 4*sin(scattering_rad*0.5)*M_PI/wavelength;
double ki = 2. * M_PI / wavelength;
dataws->dataX(iws)[0] = ki;
dataws->dataX(iws)[1] = ki + 0.00001;
Expand Down
80 changes: 80 additions & 0 deletions Code/Mantid/Framework/DataHandling/test/LoadSpiceXML2DDetTest.h
View file
Expand Up @@ -359,6 +359,86 @@ class LoadSpiceXML2DDetTest : public CxxTest::TestSuite {
AnalysisDataService::Instance().remove("Exp0335_S0038B");
AnalysisDataService::Instance().remove("Exp0335_S0038");
}


void test_LoadHB3AXMLInstrumentNoTable()
{
// Test 2theta = 15 degree
LoadSpiceXML2DDet loader;
loader.initialize();

const std::string filename("HB3A_exp355_scan0001_0522.xml");
TS_ASSERT_THROWS_NOTHING(loader.setProperty("Filename", filename));
TS_ASSERT_THROWS_NOTHING(
loader.setProperty("OutputWorkspace", "Exp0335_S0038"));
std::vector<size_t> sizelist(2);
sizelist[0] = 256;
sizelist[1] = 256;
loader.setProperty("DetectorGeometry", sizelist);
loader.setProperty("LoadInstrument", true);

loader.execute();
TS_ASSERT(loader.isExecuted());

// Get data
MatrixWorkspace_sptr outws = boost::dynamic_pointer_cast<MatrixWorkspace>(
AnalysisDataService::Instance().retrieve("Exp0335_S0038"));
TS_ASSERT(outws);
TS_ASSERT_EQUALS(outws->getNumberHistograms(), 256 * 256);

// Value
TS_ASSERT_DELTA(outws->readY(255)[0], 1.0, 0.0001);
TS_ASSERT_DELTA(outws->readY(253 * 256 + 9)[0], 1.0, 0.00001);

// Instrument
TS_ASSERT(outws->getInstrument());

Kernel::V3D source = outws->getInstrument()->getSource()->getPos();
Kernel::V3D sample = outws->getInstrument()->getSample()->getPos();
Kernel::V3D det0pos = outws->getDetector(0)->getPos();
Kernel::V3D det255pos = outws->getDetector(255)->getPos();
Kernel::V3D detlast0 = outws->getDetector(256 * 255)->getPos();
Kernel::V3D detlast = outws->getDetector(256 * 256 - 1)->getPos();
Kernel::V3D detmiddle =
outws->getDetector(256 / 2 * 256 + 256 / 2)->getPos();
Kernel::V3D detedgemiddle = outws->getDetector(256 * 256 / 2)->getPos();
std::cout << "\n";
std::cout << "Sample position: " << sample.X() << ", " << sample.Y() << ", "
<< sample.Z() << "\n";
std::cout << "Source position: " << source.X() << ", " << source.Y() << ", "
<< source.Z() << "\n";
std::cout << "Det Edge Middle: " << detedgemiddle.X() << ", "
<< detedgemiddle.Y() << ", " << detedgemiddle.Z() << "\n";
std::cout << "Det (1, 1) : " << det0pos.X() << ", " << det0pos.Y()
<< ", " << det0pos.Z() << "\n";

// center of the detector must be negative
TS_ASSERT_LESS_THAN(detmiddle.X(), 0.0);

// detector distance
double dist0 = sample.distance(det0pos);
double dist255 = sample.distance(det255pos);
double distlast = sample.distance(detlast);
double distlast0 = sample.distance(detlast0);
double distmiddle = sample.distance(detmiddle);

TS_ASSERT_DELTA(dist0, dist255, 0.0001);
TS_ASSERT_DELTA(dist0, distlast, 0.0001);
TS_ASSERT_DELTA(dist0, distlast0, 0.0001);
TS_ASSERT_DELTA(distmiddle, 0.3518, 0.000001);

// 2theta value
Kernel::V3D sample_source = sample - source;
std::cout << "Sample - Source = " << sample_source.X() << ", "
<< sample_source.Y() << ", " << sample_source.Z() << "\n";


Kernel::V3D detmid_sample = detmiddle - sample;
double twotheta_middle =
detmid_sample.angle(sample_source) * 180. / 3.1415926;
TS_ASSERT_DELTA(twotheta_middle, 42.70975, 0.02);
}

};

#endif /* MANTID_DATAHANDLING_LOADSPICEXML2DDETTEST_H_ */
27 changes: 24 additions & 3 deletions Code/Mantid/Framework/MDAlgorithms/inc/MantidMDAlgorithms/ConvertCWSDExpToMomentum.h
View file
Expand Up @@ -33,11 +33,29 @@ namespace MDAlgorithms {
File change history is stored at: <https://github.com/mantidproject/mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport ConvertCWSDExpToMomentum : API::Algorithm {
class DLLExport ConvertCWSDExpToMomentum : public API::Algorithm {
public:
ConvertCWSDExpToMomentum();
virtual ~ConvertCWSDExpToMomentum();

/// Algorithm's name
virtual const std::string name() const {
return "ConvertCWSDExpToMomentum";
}

/// Summary of algorithms purpose
virtual const std::string summary() const {
return "Load and convert a set of files in an HB3A experiment.";
}

/// Algorithm's version
virtual int version() const { return (1); }

/// Algorithm's category for identification
virtual const std::string category() const {
return "Diffraction;DataHandling\\Text";
}

private:
void init();
void exec();
Expand All @@ -48,6 +66,9 @@ class DLLExport ConvertCWSDExpToMomentum : API::Algorithm {
const detid_t &startdetid,
const int runnumber);

void convertToMomentum(const std::vector<double> &detPos, const double &wavelength,
std::vector<Mantid::coord_t> &qSample);

API::IMDEventWorkspace_sptr createExperimentMDWorkspace();

bool getInputs(std::string &errmsg);
Expand All @@ -61,17 +82,17 @@ class DLLExport ConvertCWSDExpToMomentum : API::Algorithm {
API::ITableWorkspace_sptr m_expDataTableWS;
API::ITableWorkspace_sptr m_detectorListTableWS;
API::IMDEventWorkspace_sptr m_outputWS;
Geometry::Instrument_sptr m_virtualInstrument;

Kernel::V3D m_samplePos;
Kernel::V3D m_sourcePos;

size_t m_iColFilename;
size_t m_iColStartDetID;

size_t m_numBins;

std::vector<double> m_extentMins;
std::vector<double> m_extentMaxs;
std::vector<size_t> m_numBins;
};

} // namespace MDAlgorithms
Expand Down
104 changes: 88 additions & 16 deletions Code/Mantid/Framework/MDAlgorithms/src/ConvertCWSDExpToMomentum.cpp
View file
Expand Up @@ -3,6 +3,7 @@
#include "MantidKernel/ArrayProperty.h"
#include "MantidGeometry/Instrument.h"
#include "MantidDataObjects/MDEventFactory.h"
#include "MantidAPI/ExperimentInfo.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"

using namespace Mantid::API;
Expand All @@ -13,6 +14,8 @@ using namespace Mantid::Geometry;
namespace Mantid {
namespace MDAlgorithms {

DECLARE_ALGORITHM(ConvertCWSDExpToMomentum)

//----------------------------------------------------------------------------------------------
/** Constructor
*/
Expand All @@ -36,6 +39,10 @@ void ConvertCWSDExpToMomentum::init() {
"DetectorTableWorkspace", "", Direction::Input),
"Name of table workspace containing all the detectors.");

declareProperty(new WorkspaceProperty<IMDEventWorkspace>(
"OutputWorkspace", "", Direction::Output),
"Name of MDEventWorkspace containing all experimental data.");

declareProperty(new ArrayProperty<double>("SourcePosition"),
"A vector of 3 doubles for position of source.");

Expand All @@ -44,6 +51,7 @@ void ConvertCWSDExpToMomentum::init() {

declareProperty(new ArrayProperty<double>("PixelDimension"),
"A vector of 8 doubles to determine a cubic pixel's size.");

}

/**
Expand All @@ -56,13 +64,17 @@ void ConvertCWSDExpToMomentum::exec() {
if (!inputvalid)
throw std::runtime_error(errmsg);

g_log.notice("[DB] Inputs are examined.");

// Create output MDEventWorkspace
std::vector<Kernel::V3D> vecDetPos;
std::vector<detid_t> vecDetID;
m_outputWS = createExperimentMDWorkspace();

// Add MDEventWorkspace
addMDEvents();

setProperty("OutputWorkspace", m_outputWS);

return;
}

API::IMDEventWorkspace_sptr
Expand All @@ -73,10 +85,20 @@ ConvertCWSDExpToMomentum::createExperimentMDWorkspace() {
parseDetectorTable(vec_detpos, vec_detid);

// FIXME - Should set create a solid instrument
Geometry::Instrument_sptr virtualinstrument = Geometry::ComponentHelper::createVirtualInstrument(m_sourcePos,
m_samplePos,
vec_detpos,
vec_detid);
m_virtualInstrument = Geometry::ComponentHelper::createVirtualInstrument(m_sourcePos, m_samplePos,
vec_detpos,
vec_detid);
if (!m_virtualInstrument)
throw std::runtime_error("Failed to create virtual instrument.");
g_log.notice() << "[DB] Virtual Instrument has " << m_virtualInstrument->getDetectorIDs().size()
<< "Detectors\n";


// Create workspace
// FIXME - Should I give out a better value???
m_extentMins.resize(3, -10.0);
m_extentMaxs.resize(3, 10.0);
m_numBins.resize(3, 100);

size_t m_nDimensions = 3;
IMDEventWorkspace_sptr mdws =
Expand All @@ -95,17 +117,22 @@ ConvertCWSDExpToMomentum::createExperimentMDWorkspace() {
dimensionNames[1] = "Q_sample_y";
dimensionNames[2] = "Q_sample_z";
Mantid::Kernel::SpecialCoordinateSystem coordinateSystem = Mantid::Kernel::QSample;

// Add dimensions
if (m_extentMins.size() != 3 || m_extentMaxs.size() != 3 || m_numBins.size() != 3)
throw std::runtime_error("The range and number of bins are not set up correctly to create MDEventWorkspace.");
else
for (size_t d = 0; d < 3; ++d)
g_log.debug() << "Direction " << d << ", Range = " << m_extentMins[d]
<< ", " << m_extentMaxs[d] << "\n";

for (size_t i = 0; i < m_nDimensions; ++i) {
std::string id = vec_ID[i];
std::string name = vec_name[i];
std::string name = dimensionNames[i];
// FIXME - Unit of momentum?
std::string units = "???";
// int nbins = 100;

for (size_t d = 0; d < 3; ++d)
g_log.debug() << "Direction " << d << ", Range = " << m_extentMins[d]
<< ", " << m_extentMaxs[d] << "\n";
mdws->addDimension(
Geometry::MDHistoDimension_sptr(new Geometry::MDHistoDimension(
id, name, units, static_cast<coord_t>(m_extentMins[i]),
Expand Down Expand Up @@ -158,20 +185,35 @@ void ConvertCWSDExpToMomentum::addMDEvents() {
void ConvertCWSDExpToMomentum::convertSpiceMatrixToMomentumMDEvents(
MatrixWorkspace_sptr dataws, const detid_t &startdetid,
const int runnumber) {

// Creates a new instance of the MDEventInserter.
MDEventWorkspace<MDEvent<3>, 3>::sptr mdws_mdevt_3 =
boost::dynamic_pointer_cast<MDEventWorkspace<MDEvent<3>, 3> >(m_outputWS);
MDEventInserter<MDEventWorkspace<MDEvent<3>, 3>::sptr> inserter(
mdws_mdevt_3);

size_t numspec = dataws->getNumberHistograms();
for (size_t iws = 0; iws < numspec; ++iws) {
// Get detector positions and signal
double signal = dataws->readY(iws)[0];
double error = dataws->readE(iws)[0];
double wavelength = dataws->readX(iws)[0];

// Create the MDEvent
Kernel::V3D detpos = dataws->getDetector(iws)->getPos();
std::vector<double> momentum =
convertToMomentum(detpos, wavelength, goniometersetup);
std::vector<Mantid::coord_t> momentum(3);
convertToMomentum(detpos, wavelength, momentum);
detid_t detid = dataws->getDetector(iws)->getID() + startdetid;

// Create the MDEvent
m_outputWS->addExperimentInfo();
// Insert
inserter.insertMDEvent(static_cast<float>(signal), static_cast<float>(error*error),
static_cast<uint16_t>(runnumber), detid,
momentum.data());
// m_outputWS->addExperimentInfo();
}

ExperimentInfo_sptr expinfo = boost::make_shared<ExperimentInfo>();
expinfo->setInstrument(m_virtualInstrument);
m_outputWS->addExperimentInfo(expinfo);

}

/// Examine input
Expand Down Expand Up @@ -231,6 +273,36 @@ bool ConvertCWSDExpToMomentum::getInputs(std::string &errmsg) {
return (errmsg.size() > 0);
}

void ConvertCWSDExpToMomentum::convertToMomentum(const std::vector<double> &detPos, const double &wavelength, std::vector<coord_t> &qSample)
{
// TODO/FIXME : Implement ASAP!

return;
}

API::MatrixWorkspace_sptr ConvertCWSDExpToMomentum::loadSpiceData(const std::string &filename, bool &loaded, std::string &errmsg)
{
IAlgorithm_sptr loader = createChildAlgorithm("LoadSpiceXML2DDet");
loader->initialize();
loader->setProperty("Filename", filename);
std::vector<size_t> sizelist(2);
sizelist[0] = 256;
sizelist[1] = 256;
loader->setProperty("DetectorGeometry", sizelist);
loader->setProperty("LoadInstrument", true);

loader->execute();

API::MatrixWorkspace_sptr dataws = loader->getProperty("OutputWorkspace");
if (dataws)
loaded = true;
else
loaded = false;

return dataws;
}


void ConvertCWSDExpToMomentum::parseDetectorTable(
std::vector<Kernel::V3D> &vec_detpos, std::vector<detid_t> &vec_detid) {
// Set vectors' sizes
Expand Down

0 comments on commit a2e2261

Please sign in to comment.