Skip to content

Commit

Permalink
Re #9588 Minimal changes. This is probably be deprecated in the future
Browse files Browse the repository at this point in the history
  • Loading branch information
Ricardo Ferraz Leal committed Jun 3, 2014
1 parent cdc1f46 commit 87f20df
Show file tree
Hide file tree
Showing 2 changed files with 153 additions and 145 deletions.
Expand Up @@ -43,46 +43,49 @@ namespace Algorithms {
*/
class DLLExport CorrectFlightPaths: public API::Algorithm {
public:
/// Default constructor
CorrectFlightPaths();
/// Destructor
virtual ~CorrectFlightPaths() {
}
;
/// Algorithm's name for identification overriding a virtual method
virtual const std::string name() const { return "CorrectFlightPaths"; }
/// Default constructor
CorrectFlightPaths();
/// Destructor
virtual ~CorrectFlightPaths() {
}
;
/// Algorithm's name for identification overriding a virtual method
virtual const std::string name() const {
return "CorrectFlightPaths";
}
///Summary of algorithms purpose
virtual const std::string summary() const {return "Used to correct flight paths in 2D shaped detectors.";}
/// Algorithm's version for identification overriding a virtual method
virtual int version() const {
return 1;
}
/// Algorithm's category for identification overriding a virtual method
virtual const std::string category() const {
return "Inelastic;CorrectionFunctions";
}
virtual const std::string summary() const {
return "Used to correct flight paths in 2D shaped detectors.";
}
/// Algorithm's version for identification overriding a virtual method
virtual int version() const {
return 1;
}
/// Algorithm's category for identification overriding a virtual method
virtual const std::string category() const {
return "Inelastic;CorrectionFunctions";
}

private:

// Overridden Algorithm methods
void init();
void exec();
void initWorkspaces();
double getRunProperty(std::string);
double getInstrumentProperty(std::string);
double calculateTOF(double);

/// The user selected (input) workspace
API::MatrixWorkspace_const_sptr m_inputWS;
/// The output workspace, maybe the same as the input one
API::MatrixWorkspace_sptr m_outputWS;


Geometry::Instrument_const_sptr m_instrument;
Geometry::IComponent_const_sptr m_sample;

double m_l2;
double m_wavelength;

// Overridden Algorithm methods
void init();
void exec();
void initWorkspaces();
double getRunProperty(std::string);
double getInstrumentProperty(std::string);
double calculateTOF(double);

/// The user selected (input) workspace
API::MatrixWorkspace_const_sptr m_inputWS;
/// The output workspace, maybe the same as the input one
API::MatrixWorkspace_sptr m_outputWS;

Geometry::Instrument_const_sptr m_instrument;
Geometry::IComponent_const_sptr m_sample;

double m_l2;
double m_wavelength;

};

Expand Down
221 changes: 113 additions & 108 deletions Code/Mantid/Framework/Algorithms/src/CorrectFlightPaths.cpp
@@ -1,12 +1,18 @@
/*WIKI*
Corrects the flight paths of a flat detector.
Both TOF sample-detector and distance sample-detector are corrected to constant values, i.e., this algorithm make the detector spherical rather than flat.
detector_distance must exist in the <instrument>_Parameters.xml:
Corrects the flight paths of a flat detector.
So far this has only be tested on ILL IN5.
Both TOF sample-detector and distance sample-detector are corrected to constant values, i.e., this algorithm make the detector spherical rather than flat.
*WIKI*/
The sample to detector distance must be specified as l2 in the instrument parameters file.
So far this has only be tested on the ILL IN5.
==== NOTE ====
This algorithm was coded as a proof of concept. It may be deprecated in the future.
*WIKI*/
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
Expand Down Expand Up @@ -34,7 +40,7 @@ DECLARE_ALGORITHM(CorrectFlightPaths)

// Constructor
CorrectFlightPaths::CorrectFlightPaths() :
API::Algorithm() {
API::Algorithm() {
}

/** Initialisation method. Declares properties to be used in algorithm.
Expand All @@ -44,17 +50,16 @@ 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");
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");

}

Expand All @@ -63,21 +68,20 @@ void CorrectFlightPaths::init() {
*
*/
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;
// 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;

}

Expand All @@ -87,69 +91,70 @@ void CorrectFlightPaths::initWorkspaces() {
*/
void CorrectFlightPaths::exec() {

initWorkspaces();
initWorkspaces();

Geometry::ParameterMap& pmap = m_outputWS->instrumentParameters();
Geometry::ParameterMap& pmap = m_outputWS->instrumentParameters();

const size_t numberOfChannels = this->m_inputWS->blocksize();
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
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);
// 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);

}

Expand All @@ -159,40 +164,40 @@ void CorrectFlightPaths::exec() {
*
*/
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;
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]);
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
double velocity = PhysicalConstants::h
/ (PhysicalConstants::NeutronMass * m_wavelength * 1e-10); //m/s

return distance / velocity;
return distance / velocity;
}

} // namespace Algorithm
Expand Down

0 comments on commit 87f20df

Please sign in to comment.