//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/CorrectFlightPaths.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/ListValidator.h"
#include "MantidGeometry/Instrument/ComponentHelper.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include <cmath>
namespace Mantid {
namespace Algorithms {
using namespace Kernel;
using namespace API;
using namespace Geometry;
// Register the class into the algorithm factory
DECLARE_ALGORITHM(CorrectFlightPaths)
// Constructor
CorrectFlightPaths::CorrectFlightPaths() :
API::Algorithm() {
}
/** Initialisation method. Declares properties to be used in algorithm.
*
*/
void CorrectFlightPaths::init() {
//todo: add validator for TOF
auto wsValidator = boost::make_shared<CompositeValidator>();
wsValidator->add<WorkspaceUnitValidator>("TOF");
wsValidator->add<HistogramValidator>();
declareProperty(
new WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "",
Direction::Input, wsValidator), "Name of the input workspace");
declareProperty(
new WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace", "",
Direction::Output),
"Name of the output workspace, can be the same as the input");
}
/**
* Initialises input and output workspaces.
*
*/
void CorrectFlightPaths::initWorkspaces() {
// Get the workspaces
m_inputWS = this->getProperty("InputWorkspace");
m_outputWS = this->getProperty("OutputWorkspace");
m_instrument = m_inputWS->getInstrument();
m_sample = m_instrument->getSample();
// If input and output workspaces are not the same, create a new workspace for the output
if (m_outputWS != this->m_inputWS) {
m_outputWS = API::WorkspaceFactory::Instance().create(m_inputWS);
}
m_wavelength = getRunProperty("wavelength");
g_log.debug() << "Wavelength = " << m_wavelength;
m_l2 = getInstrumentProperty("l2");
g_log.debug() << " L2 = " << m_l2 << std::endl;
}
/**
* Executes the algorithm
*
*/
void CorrectFlightPaths::exec() {
initWorkspaces();
Geometry::ParameterMap& pmap = m_outputWS->instrumentParameters();
const size_t numberOfChannels = this->m_inputWS->blocksize();
// Calculate the number of spectra in this workspace
const int numberOfSpectra = static_cast<int>(this->m_inputWS->size()
/ numberOfChannels);
API::Progress prog(this, 0.0, 1.0, numberOfSpectra);
int64_t numberOfSpectra_i = static_cast<int64_t>(numberOfSpectra); // cast to make openmp happy
// Loop over the histograms (detector spectra)
PARALLEL_FOR2(m_inputWS,m_outputWS)
for (int64_t i = 0; i < numberOfSpectra_i; ++i) {
//for (int64_t i = 32000; i < 32256; ++i) {
PARALLEL_START_INTERUPT_REGION
MantidVec& xOut = m_outputWS->dataX(i);
MantidVec& yOut = m_outputWS->dataY(i);
MantidVec& eOut = m_outputWS->dataE(i);
const MantidVec& xIn = m_inputWS->readX(i);
const MantidVec& yIn = m_inputWS->readY(i);
const MantidVec& eIn = m_inputWS->readE(i);
//Copy the energy transfer axis
// TOF
// MantidVec& xOut = m_outputWS->dataX(i);
// const MantidVec& xIn = m_inputWS->readX(i);
// subract the diference in l2
IDetector_const_sptr det = m_inputWS->getDetector(i);
double thisDetL2 = det->getDistance(*m_sample);
//if (!det->isMonitor() && thisDetL2 != m_l2) {
double deltaL2 = std::abs(thisDetL2 - m_l2);
double deltaTOF = calculateTOF(deltaL2);
deltaTOF *= 1e6; //micro sec
// position - set all detector distance to constant l2
double r, theta, phi;
V3D oldPos = det->getPos();
oldPos.getSpherical(r, theta, phi);
V3D newPos;
newPos.spherical(m_l2, theta, phi);
ComponentHelper::moveComponent(*det, pmap, newPos,
ComponentHelper::Absolute);
unsigned int j = 0;
for (; j < numberOfChannels; ++j) {
xOut[j] = xIn[j] - deltaTOF;
// there's probably a better way of copying this....
yOut[j] = yIn[j];
eOut[j] = eIn[j];
}
// last bin
xOut[numberOfChannels] = xIn[numberOfChannels] + deltaTOF;
//}
prog.report("Aligning elastic line...");
PARALLEL_END_INTERUPT_REGION
} //end for i
PARALLEL_CHECK_INTERUPT_REGION
this->setProperty("OutputWorkspace", this->m_outputWS);
}
/*
* Get run property as double
* @s - input property name
*
*/
double CorrectFlightPaths::getRunProperty(std::string s) {
Mantid::Kernel::Property* prop = m_inputWS->run().getProperty(s);
double val;
if (!prop || !Strings::convert(prop->value(), val)) {
std::string mesg = "Run property " + s + "doesn't exist!";
g_log.error(mesg);
throw std::runtime_error(mesg);
}
return val;
}
/*
* Get instrument property as double
* @s - input property name
*
*/
double CorrectFlightPaths::getInstrumentProperty(std::string s) {
std::vector<std::string> prop = m_instrument->getStringParameter(s);
if (prop.empty()) {
std::string mesg = "Property <" + s + "> doesn't exist!";
g_log.error(mesg);
throw std::runtime_error(mesg);
}
g_log.debug() << "prop[0] = " << prop[0] << std::endl;
return boost::lexical_cast<double>(prop[0]);
}
/*
* Returns the neutron TOF
* @distance - Distance in meters
*/
double CorrectFlightPaths::calculateTOF(double distance) {
double velocity = PhysicalConstants::h
/ (PhysicalConstants::NeutronMass * m_wavelength * 1e-10); //m/s
return distance / velocity;
}
} // namespace Algorithm
} // namespace Mantid
syntax highlighted by Code2HTML, v. 0.9.1