LoadSpice2D.cpp

//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidDataHandling/LoadSpice2D.h"
#include "MantidDataHandling/XmlHandler.h"
#include "MantidAPI/FileProperty.h"
#include "MantidAPI/RegisterFileLoader.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/Strings.h"
#include "MantidAPI/AlgorithmFactory.h"

#include <boost/regex.hpp>
#include <boost/shared_array.hpp>
#include <Poco/Path.h>
#include <Poco/StringTokenizer.h>
#include <Poco/DOM/DOMParser.h>
#include <Poco/DOM/Document.h>
#include <Poco/DOM/Element.h>
#include <Poco/DOM/NodeList.h>
#include <Poco/DOM/Node.h>
#include <Poco/DOM/Text.h>
#include <Poco/SAX/InputSource.h>

#include <iostream>
#include <sstream>
#include <vector>
#include <string>
#include <algorithm>

//-----------------------------------------------------------------------

using Poco::XML::DOMParser;
using Poco::XML::Document;
using Poco::XML::Element;
using Poco::XML::NodeList;
using Poco::XML::Node;
using Poco::XML::Text;

namespace Mantid {
namespace DataHandling {
// Register the algorithm into the AlgorithmFactory
DECLARE_FILELOADER_ALGORITHM(LoadSpice2D)

// Parse string and convert to numeric type
template<class T>
bool from_string(T &t, const std::string &s,
		std::ios_base &(*f)(std::ios_base &)) {
	std::istringstream iss(s);
	return !(iss >> f >> t).fail();
}


/**
 * Convenience function to store a detector value into a given spectrum.
 * Note that this type of data doesn't use TOD, so that we use a single dummy
 * bin in X. Each detector is defined as a spectrum of length 1.
 * @param ws: workspace
 * @param specID: ID of the spectrum to store the value in
 * @param value: value to store [count]
 * @param error: error on the value [count]
 * @param wavelength: wavelength value [Angstrom]
 * @param dwavelength: error on the wavelength [Angstrom]
 */
void store_value(DataObjects::Workspace2D_sptr ws, int specID, double value,
		double error, double wavelength, double dwavelength) {
	MantidVec &X = ws->dataX(specID);
	MantidVec &Y = ws->dataY(specID);
	MantidVec &E = ws->dataE(specID);
	// The following is mostly to make Mantid happy by defining a histogram with
	// a single bin around the neutron wavelength
	X[0] = wavelength - dwavelength / 2.0;
	X[1] = wavelength + dwavelength / 2.0;
	Y[0] = value;
	E[0] = error;
	ws->getSpectrum(specID)->setSpectrumNo(specID);
}

/**
 * Return the confidence with with this algorithm can load the file
 * @param descriptor A descriptor for the file
 * @returns An integer specifying the confidence level. 0 indicates it will not
 * be used
 */
int LoadSpice2D::confidence(Kernel::FileDescriptor &descriptor) const {
	if (descriptor.extension().compare(".xml") != 0)
		return 0;

	std::istream &is = descriptor.data();
	int confidence(0);

	{ // start of inner scope
		Poco::XML::InputSource src(is);
		// Set up the DOM parser and parse xml file
		DOMParser pParser;
		Poco::AutoPtr<Document> pDoc;
		try {
			pDoc = pParser.parse(&src);
		} catch (...) {
			throw Kernel::Exception::FileError("Unable to parse File:",
					descriptor.filename());
		}
		// Get pointer to root element
		Element *pRootElem = pDoc->documentElement();
		if (pRootElem) {
			if (pRootElem->tagName().compare("SPICErack") == 0) {
				confidence = 80;
			}
		}
	} // end of inner scope

	return confidence;
}

/// Constructor
LoadSpice2D::LoadSpice2D() :
		m_wavelength_input(0), m_wavelength_spread_input(0),
		m_numberXPixels(0), m_numberYPixels(0),
		m_wavelength(0), m_dwavelength(0){
}

/// Destructor
LoadSpice2D::~LoadSpice2D() {
}

/// Overwrites Algorithm Init method.
void LoadSpice2D::init() {
	declareProperty(
			new API::FileProperty("Filename", "", API::FileProperty::Load,
					".xml"), "The name of the input xml file to load");
	declareProperty(
			new API::WorkspaceProperty<API::Workspace>("OutputWorkspace", "",
					Kernel::Direction::Output),
			"The name of the Output workspace");

	// Optionally, we can specify the wavelength and wavelength spread and
	// overwrite
	// the value in the data file (used when the data file is not populated)
	auto mustBePositive =
			boost::make_shared<Kernel::BoundedValidator<double>>();
	mustBePositive->setLower(0.0);
	declareProperty("Wavelength", EMPTY_DBL(), mustBePositive,
			"Optional wavelength value to use when loading the data file "
					"(Angstrom). This value will be used instead of the value "
					"found in the data file.");
	declareProperty("WavelengthSpread", 0.1, mustBePositive,
			"Optional wavelength spread value to use when loading the "
					"data file (Angstrom). This value will be used instead of "
					"the value found in the data file.");
}


/*
 * Main method.
 * Creates an XML handler. All tag values will be a map.
 * Creates and loads the workspace with the data
 *
 */
void LoadSpice2D::exec() {

	setInputPropertiesAsMemberProperties();
	std::map<std::string, std::string> metadata = m_xmlHandler.get_metadata("Detector");
	setWavelength(metadata);

	std::vector<int> data = getData("//Data/Detector" );
	double monitorCounts = 0;
	from_string<double>(monitorCounts, metadata["Counters/monitor"], std::dec);
	double countingTime = 0;
	from_string<double>(countingTime, metadata["Counters/time"], std::dec);
	createWorkspace(data, metadata["Header/Scan_Title"], monitorCounts, countingTime );

	// Add all metadata to the WS
	addMetadataAsRunProperties(metadata);

	setMetadataAsRunProperties(metadata);

	// run load instrument
	std::string instrument = metadata["Header/Instrument"];
	runLoadInstrument(instrument, m_workspace);
	runLoadMappingTable(m_workspace, m_numberXPixels, m_numberYPixels);

	// sample_detector_distances
	double detector_distance = totalDetectorDistance(metadata);
	moveDetector(detector_distance);


}

/**
 * Parse the 2 integers of the form: INT32[192,256]
 * @param dims_str : INT32[192,256]
 */
void LoadSpice2D::parseDetectorDimensions(const std::string &dims_str) {

	// Read in the detector dimensions from the Detector tag

	boost::regex b_re_sig("INT\\d+\\[(\\d+),(\\d+)\\]");
	if (boost::regex_match(dims_str, b_re_sig)) {
		boost::match_results<std::string::const_iterator> match;
		boost::regex_search(dims_str, match, b_re_sig);
		// match[0] is the full string
		Kernel::Strings::convert(match[1], m_numberXPixels);
		Kernel::Strings::convert(match[2], m_numberYPixels);
	}

}

/**
 * Adds map of the form key:value
 * as Workspace run properties
 */
void LoadSpice2D::addMetadataAsRunProperties(
		const std::map<std::string, std::string> &metadata) {

	for (auto it = metadata.begin(); it != metadata.end(); it++) {
		std::string key = it->first;
		std::replace( key.begin(), key.end(), '/', '_');
		m_workspace->mutableRun().addProperty(key, it->second, true);
	}

}

/**
 * Get the input algorithm properties and sets them as class attributes
 */
void LoadSpice2D::setInputPropertiesAsMemberProperties() {

	m_wavelength_input = getProperty("Wavelength");
	m_wavelength_spread_input = getProperty("WavelengthSpread");

	g_log.debug() << "setInputPropertiesAsMemberProperties: "<< m_wavelength_input << " , " << m_wavelength_input << std::endl;

	std::string fileName = getPropertyValue("Filename");
	// Set up the XmlHandler handler and parse xml file
	try {
		m_xmlHandler = XmlHandler(fileName);
	} catch (...) {
		throw Kernel::Exception::FileError("Unable to parse File:", fileName);
	}
}

/**
 * Gets the wavelenght and wavelength spread from the  metadata
 * and sets them as class attributes
 */
void LoadSpice2D::setWavelength(std::map<std::string, std::string> &metadata){
	// Read in wavelength and wavelength spread

	g_log.debug() << "setWavelength: "<< m_wavelength_input << " , " << m_wavelength_input << std::endl;

	if (isEmpty(m_wavelength_input)) {
		std::string s = metadata["Header/wavelength"];
		from_string<double>(m_wavelength, s , std::dec);
		s = metadata["Header/wavelength_spread"];
		from_string<double>(m_dwavelength,s, std::dec);

		g_log.debug() << "setWavelength: "<< m_wavelength << " , " << m_dwavelength << std::endl;

	} else {
		m_wavelength = m_wavelength_input;
		m_dwavelength = m_wavelength_spread_input;
	}
}
/**
 * Parses the data dimensions and stores them as member variables
 * Reads the data and returns a vector
 */
std::vector<int> LoadSpice2D::getData(const std::string &dataXpath = "//Data/Detector" ) {
	//type : INT32[192,256]
	std::map<std::string, std::string> attributes =
			m_xmlHandler.get_attributes_from_tag(dataXpath);
	parseDetectorDimensions(attributes["type"]);
	if (m_numberXPixels == 0 || m_numberYPixels == 0)
		g_log.notice() << "Could not read in the number of pixels!"
				<< std::endl;

	std::string data_str = m_xmlHandler.get_text_from_tag(dataXpath);
	//convert string data into a vector<int>
	std::vector<int> data;
	std::stringstream iss( data_str );
	int number;
	while ( iss >> number ){
	  data.push_back( number );
	}

	if (data.size() != static_cast<size_t>(m_numberXPixels * m_numberYPixels))
		throw Kernel::Exception::NotImplementedError("Inconsistent data set: There were more data pixels found than declared in the Spice XML meta-data.");
	return data;
}

/**
 * Creates workspace and loads the data along with 2 monitors!
 */
void LoadSpice2D::createWorkspace(const std::vector<int> &data, const std::string &title,
		double monitor1_counts, double monitor2_counts) {
	int nBins = 1;
	int numSpectra = m_numberXPixels * m_numberYPixels + LoadSpice2D::nMonitors;

	m_workspace = boost::dynamic_pointer_cast<
			DataObjects::Workspace2D>(
			API::WorkspaceFactory::Instance().create("Workspace2D", numSpectra,
					nBins + 1, nBins));
	m_workspace->setTitle(title);
	m_workspace->getAxis(0)->unit() = Kernel::UnitFactory::Instance().create(
			"Wavelength");
	m_workspace->setYUnit("");
	API::Workspace_sptr workspace = boost::static_pointer_cast<API::Workspace>(
			m_workspace);
	setProperty("OutputWorkspace", workspace);

	int specID = 0;
	// Store monitor counts in the beggining
	store_value(m_workspace, specID++, monitor1_counts, monitor1_counts > 0 ? sqrt(monitor1_counts) : 0.0, m_wavelength, m_dwavelength);
	store_value(m_workspace, specID++, monitor2_counts, 0.0, m_wavelength, m_dwavelength);

	// Store detector pixels
	for(auto it = data.begin(); it != data.end(); ++it) {
		double count = static_cast<double>(*it);
		// Data uncertainties, computed according to the HFIR/IGOR reduction code
		// The following is what I would suggest instead...
		// error = count > 0 ? sqrt((double)count) : 0.0;
		double error = sqrt(0.5 + fabs(count - 0.5));
		store_value(m_workspace, specID++, count, error, m_wavelength, m_dwavelength);
	}
}

/**
 * Inserts a property in the run with a new name!
 */
template<class T>
T LoadSpice2D::addRunProperty(
		std::map<std::string, std::string> &metadata,
		const std::string &oldName, const std::string &newName, const std::string &units) {

	std::string value_str = metadata[oldName];
	T value;
	from_string<T>(value, value_str, std::dec);
	m_workspace->mutableRun().addProperty(newName, value, units, true);
	return value;
}

template<class T>
void LoadSpice2D::addRunProperty(const std::string &name, const T &value, const std::string &units) {
	m_workspace->mutableRun().addProperty(name, value, units, true);
}


/**
 * Sets the beam trap as Run Property
 */
void LoadSpice2D::setBeamTrapRunProperty(std::map<std::string, std::string> &metadata){

	// Read in beam trap positions
	double trap_pos = 0;
	from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_25mm"], std::dec);

	double beam_trap_diam = 25.4;

	double highest_trap = 0;
	from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_101mm"], std::dec);

	if (trap_pos > highest_trap) {
		highest_trap = trap_pos;
		beam_trap_diam = 101.6;
	}

	from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_50mm"], std::dec);
	if (trap_pos > highest_trap) {
		highest_trap = trap_pos;
		beam_trap_diam = 50.8;
	}

	from_string<double>(trap_pos, metadata["Motor_Positions/trap_y_76mm"], std::dec);
	if (trap_pos > highest_trap) {
		beam_trap_diam = 76.2;
	}

	addRunProperty<double>("beam-trap-diameter",beam_trap_diam,"mm");
}
void LoadSpice2D::setMetadataAsRunProperties(std::map<std::string, std::string> &metadata){
	setBeamTrapRunProperty(metadata);

	// start_time
	std::map<std::string, std::string> attributes = m_xmlHandler.get_attributes_from_tag("/");
	addRunProperty<std::string>(attributes,"start_time", "start_time","");
	addRunProperty<std::string>(attributes,"start_time", "run_start","");


	//sample thickness
	addRunProperty<double>(metadata,"Header/Sample_Thickness", "sample-thickness","mm");

	addRunProperty<double>(metadata,"Header/source_aperture_size", "source-aperture-diameter","mm");
	addRunProperty<double>(metadata,"Header/sample_aperture_size", "sample-aperture-diameter","mm");
	addRunProperty<double>(metadata,"Header/source_distance", "source-sample-distance","mm");
	addRunProperty<int>(metadata,"Motor_Positions/nguides", "number-of-guides","");

	addRunProperty<double>("wavelength",m_wavelength,"Angstrom");
	addRunProperty<double>("wavelength-spread",m_dwavelength,"Angstrom");

	addRunProperty<double>(metadata,"Counters/monitor", "monitor","");
	addRunProperty<double>(metadata,"Counters/time", "timer","sec");


}


/**
 * Calculates all the detector distances and sets them as Run properties
 * @return : total_sample_detector_distance
 */
double LoadSpice2D::totalDetectorDistance(std::map<std::string, std::string> &metadata){

	// sample_detector_distances
	double sample_detector_distance = 0;
	from_string<double>(sample_detector_distance, metadata["Motor_Positions/sample_det_dist"], std::dec);
	sample_detector_distance *= 1000.0;
	addRunProperty<double>("sample_detector_distance",sample_detector_distance,"mm");
	addRunProperty<double>("sample-detector-distance",sample_detector_distance,"mm");
	double sample_detector_distance_offset = addRunProperty<double>(metadata,"Header/tank_internal_offset", "sample-detector-distance-offset","mm");
	double sample_si_window_distance = addRunProperty<double>(metadata,"Header/sample_to_flange", "sample-si-window-distance","mm");
	double total_sample_detector_distance = sample_detector_distance + sample_detector_distance_offset + sample_si_window_distance;
	addRunProperty<double>("total-sample-detector-distance",total_sample_detector_distance,"mm");
	return total_sample_detector_distance;
}

/**
 * Places the detector at the right sample_detector_distance
 */
void LoadSpice2D::moveDetector(double sample_detector_distance) {
	// Move the detector to the right position
	API::IAlgorithm_sptr mover = createChildAlgorithm(
			"MoveInstrumentComponent");

	// Finding the name of the detector object.
	std::string detID =
			m_workspace->getInstrument()->getStringParameter("detector-name")[0];

	g_log.information("Moving " + detID);
	try {
		mover->setProperty<API::MatrixWorkspace_sptr>("Workspace", m_workspace);
		mover->setProperty("ComponentName", detID);
		mover->setProperty("Z", sample_detector_distance / 1000.0);
		mover->execute();
	} catch (std::invalid_argument &e) {
		g_log.error(
				"Invalid argument to MoveInstrumentComponent Child Algorithm");
		g_log.error(e.what());
	} catch (std::runtime_error &e) {
		g_log.error(
				"Unable to successfully run MoveInstrumentComponent Child Algorithm");
		g_log.error(e.what());
	}
}

/** Run the Child Algorithm LoadInstrument (as for LoadRaw)
 * @param inst_name :: The name written in the Nexus file
 * @param localWorkspace :: The workspace to insert the instrument into
 */
void LoadSpice2D::runLoadInstrument(const std::string &inst_name,
		DataObjects::Workspace2D_sptr localWorkspace) {

	API::IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");

	// Now execute the Child Algorithm. Catch and log any error, but don't stop.
	try {
		loadInst->setPropertyValue("InstrumentName", inst_name);
		loadInst->setProperty<API::MatrixWorkspace_sptr>("Workspace",
				localWorkspace);
		loadInst->setProperty("RewriteSpectraMap", false);
		loadInst->execute();
	} catch (std::invalid_argument &) {
		g_log.information("Invalid argument to LoadInstrument Child Algorithm");
	} catch (std::runtime_error &) {
		g_log.information(
				"Unable to successfully run LoadInstrument Child Algorithm");
	}
}

/**
 * Populate spectra mapping to detector IDs
 *
 * TODO: Get the detector size information from the workspace directly
 *
 * @param localWorkspace: Workspace2D object
 * @param nxbins: number of bins in X
 * @param nybins: number of bins in Y
 */
void LoadSpice2D::runLoadMappingTable(
		DataObjects::Workspace2D_sptr localWorkspace, int nxbins, int nybins) {
	// Get the number of monitor channels
	boost::shared_ptr<const Geometry::Instrument> instrument =
			localWorkspace->getInstrument();
	std::vector<detid_t> monitors = instrument->getMonitors();
	const int nMonitors = static_cast<int>(monitors.size());

	// Number of monitors should be consistent with data file format
	if (nMonitors != LoadSpice2D::nMonitors) {
		std::stringstream error;
		error << "Geometry error for " << instrument->getName()
				<< ": Spice data format defines " << LoadSpice2D::nMonitors
				<< " monitors, " << nMonitors << " were/was found";
		throw std::runtime_error(error.str());
	}

	// Generate mapping of detector/channel IDs to spectrum ID

	// Detector/channel counter
	int icount = 0;

	// Monitor: IDs start at 1 and increment by 1
	for (int i = 0; i < nMonitors; i++) {
		localWorkspace->getSpectrum(icount)->setDetectorID(icount + 1);
		icount++;
	}

	// Detector pixels
	for (int ix = 0; ix < nxbins; ix++) {
		for (int iy = 0; iy < nybins; iy++) {
			localWorkspace->getSpectrum(icount)->setDetectorID(
					1000000 + iy * 1000 + ix);
			icount++;
		}
	}
}

/* This method throws not found error if a element is not found in the xml file
 * @param elem :: pointer to  element
 * @param name ::  element name
 * @param fileName :: xml file name
 */
void LoadSpice2D::throwException(Poco::XML::Element *elem,
		const std::string &name, const std::string &fileName) {
	if (!elem) {
		throw Kernel::Exception::NotFoundError(
				name + " element not found in Spice XML file", fileName);
	}
}
}
}