Skip to content
Snippets Groups Projects
Commit 0225af11 authored by WHITFIELDRE email's avatar WHITFIELDRE email
Browse files

Merge pull request #13126 from mantidproject/13087_ModeratorTzero_Extend 

13087 moderator tzero extend
parents a9c09da4 7397640b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ModeratorTzero.h
+ 9
− 9
View file @ 0225af11
......@@ -104,15 +104,15 @@ private:
/// Execution code for histogram workspace
void exec();
/// Execution code for event workspace
void execEvent();
/// Calculate distance from source to sample or monitor
double CalculateL1(Mantid::API::MatrixWorkspace_sptr inputWS, size_t i);
/// Calculate time from sample to detector
double CalculateT2(Mantid::API::MatrixWorkspace_sptr inputWS, size_t i);
/// Calculate emission time from the moderator for a given detector (L1, t2)
/// and TOF
double CalculateT0(const double &tof, const double &L1, const double &t2,
double &E1, mu::Parser &parser);
void execEvent(const std::string &emode);
/// Calculate emission time from the moderator for a given
/// detector (L1, t2) and TOF when Emode==Inelastic
double CalculateT0indirect(const double &tof, const double &L1,
const double &t2, double &E1, mu::Parser &parser);
/// Calculate emission time from the moderator for a given
/// detector (L1, t2) and TOF when Emode==Elastic
double CalculateT0elastic(const double &tof, const double &L12,
double &E1, mu::Parser &parser);
const double m_convfactor;
/// Maximum number of iterations when calculating the emission time from the
/// moderator
......
Code/Mantid/Framework/Algorithms/src/ModeratorTzero.cpp
+ 367
− 177
View file @ 0225af11
......@@ -5,6 +5,7 @@
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidDataObjects/Workspace2D.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidKernel/ListValidator.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidGeometry/muParser_Silent.h"
#include <boost/lexical_cast.hpp>
......@@ -37,21 +38,29 @@ void ModeratorTzero::init() {
auto wsValidator = boost::make_shared<CompositeValidator>();
wsValidator->add<WorkspaceUnitValidator>("TOF");
declareProperty(new WorkspaceProperty<MatrixWorkspace>(
"InputWorkspace", "", Direction::Input, wsValidator),
"The name of the input workspace, containing events and/or "
"histogram data, in units of time-of-flight");
"InputWorkspace", "", Direction::Input, wsValidator),
"The name of the input workspace, containing events and/or "
"histogram data, in units of time-of-flight");
// declare the output workspace
declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace", "",
Direction::Output),
"The name of the output workspace");
Direction::Output), "The name of the output workspace");
// declare the instrument scattering mode
std::vector<std::string> EModeOptions;
EModeOptions.push_back("Indirect");
EModeOptions.push_back("Direct");
EModeOptions.push_back("Elastic");
this->declareProperty("EMode", "Indirect",
boost::make_shared<StringListValidator>(EModeOptions),
"The energy mode (default: Indirect)");
declareProperty(new Kernel::PropertyWithValue<double>(
"tolTOF", 0.1, Kernel::Direction::Input),
"Tolerance in the calculation of the emission time, in "
"microseconds (default:1)");
"tolTOF", 0.1, Kernel::Direction::Input),
"Tolerance in the calculation of the emission time, in "
"microseconds (default:1)");
declareProperty(new Kernel::PropertyWithValue<size_t>(
"Niter", 1, Kernel::Direction::Input),
"Number of iterations (default:1)");
"Niter", 1, Kernel::Direction::Input),
"Number of iterations (default:1)");
} // end of void ModeratorTzero::init()
void ModeratorTzero::exec() {
......@@ -60,34 +69,52 @@ void ModeratorTzero::exec() {
m_niter = getProperty("Niter"); // number of iterations
const MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
m_instrument = inputWS->getInstrument(); // pointer to the instrument
const std::string emode = getProperty("EMode");
// deltaE-mode (should be "indirect")
std::vector<std::string> Emode =
m_instrument->getStringParameter("deltaE-mode");
if (Emode.empty())
throw Exception::InstrumentDefinitionError(
"Unable to retrieve instrument geometry (direct or indirect) parameter",
inputWS->getTitle());
if (Emode[0] != "indirect")
throw Exception::InstrumentDefinitionError(
"Instrument geometry must be of type indirect.");
// Check if Ei stored in workspace
if(emode=="Direct")
{
if(!inputWS->run().hasProperty("Ei"))
{
g_log.error("No incident energy value (Ei) has been set or stored.");
return;
}
}
// extract formula from instrument parameters
std::vector<std::string> t0_formula =
m_instrument->getStringParameter("t0_formula");
m_instrument->getStringParameter("t0_formula");
if (t0_formula.empty())
throw Exception::InstrumentDefinitionError(
"Unable to retrieve t0_formula among instrument parameters");
{
g_log.error("Unable to retrieve t0_formula among instrument parameters.");
return;
}
m_formula = t0_formula[0];
// Are there source and sample?
IComponent_const_sptr source;
IComponent_const_sptr sample;
double Lss(0); // distance from source to sample
try
{
source = m_instrument->getSource();
sample = m_instrument->getSample();
Lss = source->getDistance(*sample);
}
catch(Exception::NotFoundError &)
{
g_log.error("Unable to calculate source-sample distance.");
return;
}
// Run execEvent if eventWorkSpace
EventWorkspace_const_sptr eventWS =
boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
if (eventWS != NULL) {
execEvent();
execEvent(emode);
return;
}
// Run exec for matrix workspace
MatrixWorkspace_sptr outputWS = getProperty("OutputWorkspace");
// Check whether input == output to see whether a new workspace is required.
if (outputWS != inputWS) {
......@@ -95,81 +122,173 @@ void ModeratorTzero::exec() {
outputWS = WorkspaceFactory::Instance().create(inputWS);
}
// calculate tof shift once for all neutrons if emode==Direct
double t0_direct(-1);
if(emode=="Direct")
{
Kernel::Property *eiprop = inputWS->run().getProperty("Ei");
double Ei = boost::lexical_cast<double>(eiprop->value());
mu::Parser parser;
parser.DefineVar("incidentEnergy", &Ei); // associate E1 to this parser
parser.SetExpr(m_formula);
t0_direct = parser.Eval();
}
const size_t numHists = static_cast<size_t>(inputWS->getNumberHistograms());
Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
PARALLEL_FOR2(inputWS, outputWS)
// iterate over the spectra
for (int i = 0; i < static_cast<int>(numHists); ++i) {
for (int i = 0; i < static_cast<int>(numHists); ++i)
{
PARALLEL_START_INTERUPT_REGION
size_t wsIndex = static_cast<size_t>(i);
double L1 = CalculateL1(
inputWS, wsIndex); // distance from source to sample or monitor
double t2 = CalculateT2(inputWS, wsIndex); // time from sample to detector
// shift the time of flights by the emission time from the moderator
if (t2 >= 0) // t2 < 0 when no detector info is available
MantidVec &inbins = inputWS->dataX(i);
MantidVec &outbins = outputWS->dataX(i);
// One parser for each parallel processor needed (except Edirect mode)
double E1;
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
parser.SetExpr(m_formula);
IDetector_const_sptr det;
double L1(Lss); // distance from source to sample
double L2(-1); // distance from sample to detector
try
{
double E1;
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
parser.SetExpr(m_formula);
E1 = m_convfactor * (L1 / m_t1min) * (L1 / m_t1min);
double min_t0_next = parser.Eval(); // fast neutrons are shifted by
// min_t0_next, irrespective of tof
MantidVec &inbins = inputWS->dataX(i);
MantidVec &outbins = outputWS->dataX(i);
// iterate over the time-of-flight values
for (unsigned int ibin = 0; ibin < inbins.size(); ibin++) {
double tof = inbins[ibin]; // current time-of-flight
if (tof < m_t1min + t2)
tof -= min_t0_next;
else
tof -= CalculateT0(tof, L1, t2, E1, parser);
outbins[ibin] = tof;
det = inputWS->getDetector(i);
if (det->isMonitor())
{
// redefine the sample as the monitor
L1 = source->getDistance(*det);
L2 = 0;
}
} else {
outputWS->dataX(i) = inputWS->dataX(i);
else
{
L2 = sample->getDistance(*det);
}
} // end of try
catch(Exception::NotFoundError &)
{
g_log.error() << "Unable to calculate distances to/from detector" << i << std::endl;
outbins = inbins;
}
if(L2 >= 0)
{
// fast neutrons are shifted by min_t0_next, irrespective of tof
double v1_max = L1 / m_t1min;
E1 = m_convfactor * v1_max*v1_max;
double min_t0_next = parser.Eval();
if(emode=="Indirect")
{
double t2(-1.0); // time from sample to detector. (-1) signals error
if(det->isMonitor())
{
t2=0.0;
}
else
{
static const double convFact = 1.0e-6 * sqrt(2 * PhysicalConstants::meV /
PhysicalConstants::NeutronMass);
std::vector<double> wsProp = det->getNumberParameter("Efixed");
if (!wsProp.empty())
{
double E2 = wsProp.at(0); //[E2]=meV
double v2 = convFact * sqrt(E2); //[v2]=meter/microsec
t2 = L2 / v2;
}
else
{
// t2 is kept to -1 if no Efixed is found
g_log.debug() << "Efixed not found for detector " << i << std::endl;
}
}
// shift the time of flights by the emission time from the moderator
if (t2 >= 0) // t2 < 0 when no detector info is available
{
// iterate over the time-of-flight values
for (unsigned int ibin = 0; ibin < inbins.size(); ibin++) {
double tof = inbins[ibin]; // recorded time-of-flight
if (tof < m_t1min + t2)
tof -= min_t0_next;
else
tof -= CalculateT0indirect(tof, L1, t2, E1, parser);
outbins[ibin] = tof;
}
}
else
{
outbins = inbins;
}
} // end of if(emode=="Indirect")
else if(emode=="Elastic")
{
for (unsigned int ibin = 0; ibin < inbins.size(); ibin++)
{
double tof = inbins[ibin]; // recorded time-of-flight;
if (tof < m_t1min *(L1+L2)/L1)
tof -= min_t0_next;
else
tof -= CalculateT0elastic(tof, L1+L2, E1, parser);
outbins[ibin] = tof;
}
} // end of else if(emode=="Elastic")
else if(emode=="Direct")
{
for (unsigned int ibin = 0; ibin < inbins.size(); ibin++)
{
outbins[ibin] = inbins[ibin] - t0_direct;
}
} // end of else if(emode="Direct")
} // end of if(L2 >= 0)
// Copy y and e data
outputWS->dataY(i) = inputWS->dataY(i);
outputWS->dataE(i) = inputWS->dataE(i);
prog.report();
PARALLEL_END_INTERUPT_REGION
}
} // end of for (int i = 0; i < static_cast<int>(numHists); ++i)
PARALLEL_CHECK_INTERUPT_REGION
// Copy units
if (inputWS->getAxis(0)->unit().get()) {
if (inputWS->getAxis(0)->unit().get())
{
outputWS->getAxis(0)->unit() = inputWS->getAxis(0)->unit();
}
try {
if (inputWS->getAxis(1)->unit().get()) {
try
{
if (inputWS->getAxis(1)->unit().get())
{
outputWS->getAxis(1)->unit() = inputWS->getAxis(1)->unit();
}
} catch (Exception::IndexError &) {
}
catch (Exception::IndexError &)
{
// OK, so this isn't a Workspace2D
}
// Assign it to the output workspace property
setProperty("OutputWorkspace", outputWS);
}
} //end of void ModeratorTzero::exec
void ModeratorTzero::execEvent() {
void ModeratorTzero::execEvent(const std::string &emode)
{
g_log.information("Processing event workspace");
const MatrixWorkspace_const_sptr matrixInputWS =
getProperty("InputWorkspace");
const MatrixWorkspace_const_sptr matrixInputWS = getProperty("InputWorkspace");
EventWorkspace_const_sptr inputWS =
boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);
boost::dynamic_pointer_cast<const EventWorkspace>(matrixInputWS);
// generate the output workspace pointer
const size_t numHists = static_cast<size_t>(inputWS->getNumberHistograms());
Mantid::API::MatrixWorkspace_sptr matrixOutputWS =
getProperty("OutputWorkspace");
Mantid::API::MatrixWorkspace_sptr matrixOutputWS = getProperty("OutputWorkspace");
EventWorkspace_sptr outputWS;
if (matrixOutputWS == matrixInputWS) {
if (matrixOutputWS == matrixInputWS)
{
outputWS = boost::dynamic_pointer_cast<EventWorkspace>(matrixOutputWS);
} else {
}
else
{
// Make a brand new EventWorkspace
outputWS = boost::dynamic_pointer_cast<EventWorkspace>(
WorkspaceFactory::Instance().create("EventWorkspace", numHists, 2, 1));
......@@ -182,143 +301,194 @@ void ModeratorTzero::execEvent() {
setProperty("OutputWorkspace", matrixOutputWS);
}
// Get a pointer to the sample
IComponent_const_sptr sample = outputWS->getInstrument()->getSample();
// Get pointers to sample and source
IComponent_const_sptr source = m_instrument->getSource();
IComponent_const_sptr sample = m_instrument->getSample();
double Lss = source->getDistance(*sample); // distance from source to sample
// calculate tof shift once for all neutrons if emode==Direct
double t0_direct(-1);
if(emode=="Direct")
{
Kernel::Property *eiprop = inputWS->run().getProperty("Ei");
double Ei = boost::lexical_cast<double>(eiprop->value());
mu::Parser parser;
parser.DefineVar("incidentEnergy", &Ei); // associate E1 to this parser
parser.SetExpr(m_formula);
t0_direct = parser.Eval();
}
// Loop over the spectra
Progress prog(this, 0.0, 1.0, numHists); // report progress of algorithm
PARALLEL_FOR1(outputWS)
for (int i = 0; i < static_cast<int>(numHists); ++i) {
for (int i = 0; i < static_cast<int>(numHists); ++i)
{
PARALLEL_START_INTERUPT_REGION
size_t wsIndex = static_cast<size_t>(i);
EventList &evlist = outputWS->getEventList(wsIndex);
if (evlist.getNumberEvents() > 0) // don't bother with empty lists
{
double L1 = CalculateL1(
matrixOutputWS, wsIndex); // distance from source to sample or monitor
double t2 =
CalculateT2(matrixOutputWS, wsIndex); // time from sample to detector
if (t2 >= 0) // t2 < 0 when no detector info is available
IDetector_const_sptr det;
double L1(Lss); // distance from source to sample
double L2(-1); // distance from sample to detector
try
{
det = inputWS->getDetector(i);
if (det->isMonitor())
{
// redefine the sample as the monitor
L1 = source->getDistance(*det);
L2 = 0;
}
else
{
L2 = sample->getDistance(*det);
}
}
catch(Exception::NotFoundError &)
{
g_log.error() << "Unable to calculate distances to/from detector" << i << std::endl;
}
if(L2 >= 0)
{
double tof, E1;
// One parser for each parallel processor needed (except Edirect mode)
double E1;
mu::Parser parser;
parser.DefineVar("incidentEnergy",
&E1); // associate variable E1 to this parser
parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
parser.SetExpr(m_formula);
E1 = m_convfactor * (L1 / m_t1min) * (L1 / m_t1min);
double min_t0_next = parser.Eval(); // fast neutrons are shifted by
// min_t0_next, irrespective of tof
// fix the histogram bins
MantidVec &x = evlist.dataX();
for (MantidVec::iterator iter = x.begin(); iter != x.end(); ++iter) {
tof = *iter;
if (tof < m_t1min + t2)
tof -= min_t0_next;
else
tof -= CalculateT0(tof, L1, t2, E1, parser);
*iter = tof;
}
MantidVec tofs = evlist.getTofs();
for (unsigned int itof = 0; itof < tofs.size(); itof++) {
tof = tofs[itof] + 0.002 * (rand() % 100 - 50); // add a [-0.1,0.1]
// microsecond noise
// to avoid artifacts
// resulting from
// original tof data
if (tof < m_t1min + t2)
tof -= min_t0_next;
// fast neutrons are shifted by min_t0_next, irrespective of tof
double v1_max = L1 / m_t1min;
E1 = m_convfactor * v1_max*v1_max;
double min_t0_next = parser.Eval();
if(emode=="Indirect")
{
double t2(-1.0); // time from sample to detector. (-1) signals error
if(det->isMonitor())
{
t2=0.0;
}
else
tof -= CalculateT0(tof, L1, t2, E1, parser);
tofs[itof] = tof;
}
evlist.setTofs(tofs);
evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
}
}
{
static const double convFact = 1.0e-6 * sqrt(2 * PhysicalConstants::meV /
PhysicalConstants::NeutronMass);
std::vector<double> wsProp = det->getNumberParameter("Efixed");
if (!wsProp.empty())
{
double E2 = wsProp.at(0); //[E2]=meV
double v2 = convFact * sqrt(E2); //[v2]=meter/microsec
t2 = L2 / v2;
}
else
{
// t2 is kept to -1 if no Efixed is found
g_log.debug() << "Efixed not found for detector " << i << std::endl;
}
}
if(t2 >= 0) // t2 < 0 when no detector info is available
{
double tof;
// fix the histogram bins
MantidVec &x = evlist.dataX();
for (MantidVec::iterator iter = x.begin(); iter != x.end(); ++iter)
{
tof = *iter;
if (tof < m_t1min + t2)
tof -= min_t0_next;
else
tof -= CalculateT0indirect(tof, L1, t2, E1, parser);
*iter = tof;
}
MantidVec tofs = evlist.getTofs();
for (unsigned int itof = 0; itof < tofs.size(); itof++)
{
tof = tofs[itof];
if (tof < m_t1min + t2)
tof -= min_t0_next;
else
tof -= CalculateT0indirect(tof, L1, t2, E1, parser);
tofs[itof] = tof;
}
evlist.setTofs(tofs);
evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
} // end of if( t2>= 0)
}// end of if(emode=="Indirect")
else if(emode=="Elastic")
{
double tof;
// Apply t0 correction to histogram bins
MantidVec &x = evlist.dataX();
for (MantidVec::iterator iter = x.begin(); iter != x.end(); ++iter)
{
tof = *iter;
if(tof < m_t1min *(L1+L2)/L1)
tof -= min_t0_next;
else
tof -= CalculateT0elastic(tof, L1+L2, E1, parser);
*iter = tof;
}
MantidVec tofs = evlist.getTofs();
for (unsigned int itof = 0; itof < tofs.size(); itof++)
{
// add a [-0.1,0.1] microsecond noise to avoid artifacts
// resulting from original tof data
tof = tofs[itof];
if(tof < m_t1min *(L1+L2)/L1)
tof -= min_t0_next;
else
tof -= CalculateT0elastic(tof, L1+L2, E1, parser);
tofs[itof] = tof;
}
evlist.setTofs(tofs);
evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
MantidVec tofs_b=evlist.getTofs();
MantidVec xarray=evlist.readX();
} // end of else if(emode=="Elastic")
else if(emode=="Direct")
{
// fix the histogram bins
MantidVec &x = evlist.dataX();
for (MantidVec::iterator iter = x.begin(); iter != x.end(); ++iter)
{
*iter -= t0_direct;
}
MantidVec tofs = evlist.getTofs();
for (unsigned int itof = 0; itof < tofs.size(); itof++)
{
tofs[itof] -=t0_direct;
}
evlist.setTofs(tofs);
evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
} // end of else if(emode=="Direct")
} // end of if(L2 >= 0)
} // end of if (evlist.getNumberEvents() > 0)
prog.report();
PARALLEL_END_INTERUPT_REGION
}
} // end of for (int i = 0; i < static_cast<int>(numHists); ++i)
PARALLEL_CHECK_INTERUPT_REGION
outputWS->clearMRU(); // Clears the Most Recent Used lists */
} // end of void ModeratorTzero::execEvent()
/// calculate distance from source to sample or detector
double ModeratorTzero::CalculateL1(Mantid::API::MatrixWorkspace_sptr inputWS,
size_t i) {
double L1(0);
// Get detector position
IDetector_const_sptr det;
try {
det = inputWS->getDetector(i);
} catch (Exception::NotFoundError &) {
return 0;
}
if (det->isMonitor()) {
L1 = m_instrument->getSource()->getDistance(*det);
} else {
IComponent_const_sptr sample = m_instrument->getSample();
try {
L1 = m_instrument->getSource()->getDistance(*sample);
} catch (Exception::NotFoundError &) {
g_log.error("Unable to calculate source-sample distance");
throw Exception::InstrumentDefinitionError(
"Unable to calculate source-sample distance", inputWS->getTitle());
}
}
return L1;
}
// calculate time from sample to detector
double ModeratorTzero::CalculateT2(MatrixWorkspace_sptr inputWS, size_t i) {
static const double convFact = 1.0e-6 * sqrt(2 * PhysicalConstants::meV /
PhysicalConstants::NeutronMass);
double t2(-1.0); // negative initialization signals error
// Get detector position
IDetector_const_sptr det;
try {
det = inputWS->getDetector(i);
} catch (Exception::NotFoundError &) {
return t2;
}
/// Calculate emission time for a given detector (L1, t2)
/// and TOF when Emode==Inelastic
double ModeratorTzero::CalculateT0indirect(const double &tof, const double &L1,
const double &t2, double &E1, mu::Parser &parser) {
if (det->isMonitor()) {
t2 = 0.0; // t2=0.0 since there is no sample to detector path
} else {
IComponent_const_sptr sample = m_instrument->getSample();
// Get final energy E_f, final velocity v_f
std::vector<double> wsProp = det->getNumberParameter("Efixed");
if (!wsProp.empty()) {
double E2 = wsProp.at(0); //[E2]=meV
double v2 = convFact * sqrt(E2); //[v2]=meter/microsec
try {
double L2 = det->getDistance(*sample);
t2 = L2 / v2;
} catch (Exception::NotFoundError &) {
g_log.error("Unable to calculate detector-sample distance");
throw Exception::InstrumentDefinitionError(
"Unable to calculate detector-sample distance",
inputWS->getTitle());
}
} else {
g_log.debug() << "Efixed not found for detector " << i << std::endl;
}
}
return t2;
} // end of CalculateT2(const MatrixWorkspace_sptr inputWS, size_t i)
/// Calculate emission time for a given detector (L1, t2) and TOF
double ModeratorTzero::CalculateT0(const double &tof, const double &L1,
const double &t2, double &E1,
mu::Parser &parser) {
double t0_curr, t0_next;
t0_curr = m_tolTOF; // current iteration emission time
t0_next = 0.0; // next iteration emission time, initialized to zero
size_t iiter(0); // current iteration number
// iterate until convergence in t0 reached
while (std::fabs(t0_curr - t0_next) >= m_tolTOF && iiter < m_niter) {
while (std::fabs(t0_curr - t0_next) >= m_tolTOF && iiter < m_niter)
{
t0_curr = t0_next;
double t1 = tof - t0_curr - t2;
double v1 = L1 / t1;
......@@ -329,6 +499,26 @@ double ModeratorTzero::CalculateT0(const double &tof, const double &L1,
return t0_next;
}
/// Calculate emission time for a given detector (L1, t2)
/// and TOF when Emode==Elastic
double ModeratorTzero::CalculateT0elastic(const double &tof, const double &L12,
double &E1, mu::Parser &parser)
{
double t0_curr, t0_next;
t0_curr = m_tolTOF; // current iteration emission time
t0_next = 0.0; // next iteration emission time, initialized to zero
size_t iiter(0); // current iteration number
// iterate until convergence in t0 reached
while (std::fabs(t0_curr - t0_next) >= m_tolTOF && iiter < m_niter)
{
t0_curr = t0_next;
double v1 = L12 / (tof - t0_curr); // v1 = v2 = v since emode is elastic
E1 = m_convfactor * v1 * v1; // Energy in meV if v1 in meter/microsecond
t0_next = parser.Eval();
iiter++;
}
return t0_next;
}
double ModeratorTzero::gett1min() { return m_t1min; }
} // namespace Algorithms
......
Code/Mantid/Framework/Algorithms/test/ModeratorTzeroTest.h
+ 156
− 39
View file @ 0225af11
......@@ -37,92 +37,197 @@ void TestInit()
TS_ASSERT( alg.isInitialized() );
}
void TestExecThrowsDeltaEmode()
void TestExecHistogramIndirect()
{
MatrixWorkspace_sptr testWS = CreateHistogramWorkspace();
AnalysisDataService::Instance().add("testWS", testWS);
const bool add_Efixed=true;
const bool add_t0_formula=true;
AddToIndirectInstrument(testWS, add_t0_formula, add_Efixed);
ModeratorTzero alg;
alg.initialize();
alg.setProperty("InputWorkspace",testWS);
alg.setProperty("OutputWorkspace","testWS");
alg.setRethrows(true); // necessary, otherwise the algorithm will catch all exceptions and not return them
TS_ASSERT_THROWS(alg.execute(), Exception::InstrumentDefinitionError);
alg.setProperty("EMode", "Indirect");
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check a few values. These are values separated by 400 bins
const size_t jump(400);
double tofs[3][11]={
{-0.218694, 1599.78, 3199.78, 4799.78, 6399.78, 7999.78, 9550.71, 11150.2, 12750.1, 14350, 15950},
{-34.9412, 1550.24, 3150.06, 4750.03, 6350.01, 7950.01, 9550.01, 11150, 12750, 14350, 15950},
{-9.67714, 1550.63, 3150.16, 4750.07, 6350.04, 7950.03, 9550.02, 11150, 12750, 14350, 15950}
};
for (size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
MantidVec xarray=testWS->dataX(ihist);
for(size_t ibin=0; ibin < xarray.size(); ibin+=jump)
{
TS_ASSERT_DELTA(tofs[ihist][ibin/jump],xarray[ibin],0.1);
}
}
AnalysisDataService::Instance().remove("testWS");
}
void TestExecThrowsNoFormula()
void TestExecHistogramElastic()
{
MatrixWorkspace_sptr testWS = CreateHistogramWorkspace();
AnalysisDataService::Instance().add("testWS", testWS);
const bool add_deltaE_mode=true;
AddToInstrument(testWS,add_deltaE_mode);
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),
"t0_formula","101.9*incidentEnergy^(-0.41)*exp(-incidentEnergy/282.0)");
ModeratorTzero alg;
alg.initialize();
alg.setProperty("InputWorkspace",testWS);
alg.setProperty("OutputWorkspace","testWS");
alg.setRethrows(true); // necessary, otherwise the algorithm will catch all exceptions and not return them
TS_ASSERT_THROWS(alg.execute(), Exception::InstrumentDefinitionError);
alg.setProperty("EMode", "Elastic");
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check a few values. These are values separated by 400 bins
const size_t jump(400);
double tofs[3][11]={
{0.0, 1599.94, 3196.91, 4791.92, 6387.25, 7983.04, 9579.19, 11175.6, 12772.2, 14368.9, 15965.7},
{0.0, 1595.57, 3184.9, 4776.23, 6368.59, 7961.54, 9554.85, 11148.4, 12742.2, 14336.2, 15930.3},
{0.0, 1599.32, 3193.02, 4786.52, 6380.87, 7975.78, 9571.06, 11166.6, 12762.3, 14358.2, 15954.1}
};
for (size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
MantidVec xarray=testWS->dataX(ihist);
for(size_t ibin=0; ibin < xarray.size(); ibin+=jump)
{
TS_ASSERT_DELTA(tofs[ihist][ibin/jump],xarray[ibin],0.1);
}
}
AnalysisDataService::Instance().remove("testWS");
}
/*
* First spectrum is a detector. Remaining two spectra are monitors
*/
void TestExecHistogram()
void TestExecHistogramDirect()
{
MatrixWorkspace_sptr testWS = CreateHistogramWorkspace();
AnalysisDataService::Instance().add("testWS", testWS);
const bool add_deltaE_mode=true;
const bool add_t0_formula=true;
AddToInstrument(testWS, add_deltaE_mode, add_t0_formula);
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),
"t0_formula","101.9*incidentEnergy^(-0.41)*exp(-incidentEnergy/282.0)");
testWS->mutableRun().addProperty("Ei", 100.0, "meV", true);
ModeratorTzero alg;
alg.initialize();
alg.setProperty("InputWorkspace",testWS);
alg.setProperty("OutputWorkspace","testWS");
alg.setProperty("EMode", "Direct");
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check a few values
// Check a few values. These are values separated by 400 bins
const size_t jump(400);
double tofs[3][11]={
{-0.218694, 1599.78, 3199.78, 4799.78, 6399.78, 7999.78, 9550.71, 11150.2, 12750.1, 14350, 15950},
{-34.9412, 1550.24, 3150.06, 4750.03, 6350.01, 7950.01, 9550.01, 11150, 12750, 14350, 15950},
{-9.67714, 1550.63, 3150.16, 4750.07, 6350.04, 7950.03, 9550.02, 11150, 12750, 14350, 15950}
{-10.8185,1589.18,3189.18,4789.18,6389.18,7989.18,9589.18,11189.2,12789.2,14389.2,15989.2},
{-10.8185,1589.18,3189.18,4789.18,6389.18,7989.18,9589.18,11189.2,12789.2,14389.2,15989.2},
{-10.8185,1589.18,3189.18,4789.18,6389.18,7989.18,9589.18,11189.2,12789.2,14389.2,15989.2}
};
for (size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
MantidVec xarray=testWS->dataX(ihist);
for(size_t ibin=0; ibin < xarray.size(); ibin+=400)
TS_ASSERT_DELTA(tofs[ihist][ibin/400],xarray[ibin],0.1);
for(size_t ibin=0; ibin < xarray.size(); ibin+=jump)
{
TS_ASSERT_DELTA(tofs[ihist][ibin/jump],xarray[ibin],0.1);
}
}
AnalysisDataService::Instance().remove("testWS");
}
void TestExecEvents()
void TestExecEventsIndirect()
{
EventWorkspace_sptr testWS=CreateEventWorkspace();
AnalysisDataService::Instance().add("testWS", testWS);
const bool add_deltaE_mode=true;
const bool add_Efixed=true;
const bool add_t0_formula=true;
MatrixWorkspace_sptr mtestWS=boost::dynamic_pointer_cast<MatrixWorkspace>(testWS);
AddToInstrument(mtestWS, add_deltaE_mode, add_t0_formula);
AddToIndirectInstrument(mtestWS, add_t0_formula, add_Efixed);
ModeratorTzero alg;
alg.initialize();
alg.setProperty("InputWorkspace",testWS);
alg.setProperty("OutputWorkspace","testWS");
alg.setProperty("EMode", "Indirect");
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check a few values
// Check a few values. These are values separated by 400 bins
const size_t jump(400);
double tofs_a[11]={-37.5547, 1562.45, 3162.45, 4762.45, 6362.45, 7962.45, 9550.18, 11150, 12750, 14350, 15950};
for (size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
EventList &evlist=testWS->getEventList(ihist);
MantidVec tofs_b=evlist.getTofs();
MantidVec xarray=evlist.readX();
for(size_t ibin=0; ibin < xarray.size(); ibin+=400)
for(size_t ibin=0; ibin < xarray.size(); ibin+=jump)
{
TS_ASSERT_DELTA(tofs_a[ibin/400],xarray[ibin],0.1);
TS_ASSERT_DELTA(tofs_a[ibin/400],tofs_b[ibin],0.2);
TS_ASSERT_DELTA(tofs_a[ibin/jump],xarray[ibin],0.1);
TS_ASSERT_DELTA(tofs_a[ibin/jump],tofs_b[ibin],0.2);
}
}
AnalysisDataService::Instance().remove("testWS");
}
void TestExecEventsElastic()
{
EventWorkspace_sptr testWS=CreateEventWorkspace();
AnalysisDataService::Instance().add("testWS", testWS);
MatrixWorkspace_sptr mtestWS=boost::dynamic_pointer_cast<MatrixWorkspace>(testWS);
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),
"t0_formula","101.9*incidentEnergy^(-0.41)*exp(-incidentEnergy/282.0)");
ModeratorTzero alg;
alg.initialize();
alg.setProperty("InputWorkspace",testWS);
alg.setProperty("OutputWorkspace","testWS");
alg.setProperty("EMode", "Elastic");
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check a few values. These are values separated by 400 bins
const size_t jump(400);
double tofs_a[11]={0.0,1598.38,3190.3,4783.04,6376.76,7971.06,9565.72,11160.6,12755.7,14351,15946.3};
for (size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
EventList &evlist=testWS->getEventList(ihist);
MantidVec tofs_b=evlist.getTofs();
MantidVec xarray=evlist.readX();
for(size_t ibin=0; ibin < xarray.size(); ibin+=jump)
{
TS_ASSERT_DELTA(tofs_a[ibin/jump],xarray[ibin],0.1);
TS_ASSERT_DELTA(tofs_a[ibin/jump],tofs_b[ibin],0.2);
}
}
AnalysisDataService::Instance().remove("testWS");
}
void TestExecEventsDirect()
{
EventWorkspace_sptr testWS=CreateEventWorkspace();
AnalysisDataService::Instance().add("testWS", testWS);
MatrixWorkspace_sptr mtestWS=boost::dynamic_pointer_cast<MatrixWorkspace>(testWS);
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),
"t0_formula","101.9*incidentEnergy^(-0.41)*exp(-incidentEnergy/282.0)");
testWS->mutableRun().addProperty("Ei", 100.0, "meV", true);
ModeratorTzero alg;
alg.initialize();
alg.setProperty("InputWorkspace",testWS);
alg.setProperty("OutputWorkspace","testWS");
alg.setProperty("EMode", "Direct");
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check a few values. These are values separated by 400 bins
const size_t jump(400);
double tofs_a[11]={-10.8185,1589.18,3189.18,4789.18,6389.18,7989.18,9589.18,11189.2,12789.2,14389.2,15989.2};
for (size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
EventList &evlist=testWS->getEventList(ihist);
MantidVec tofs_b=evlist.getTofs();
MantidVec xarray=evlist.readX();
for(size_t ibin=0; ibin < xarray.size(); ibin+=jump)
{
TS_ASSERT_DELTA(tofs_a[ibin/jump],xarray[ibin],0.1);
TS_ASSERT_DELTA(tofs_a[ibin/jump],tofs_b[ibin],0.2);
}
}
AnalysisDataService::Instance().remove("testWS");
......@@ -130,6 +235,10 @@ void TestExecEvents()
private:
/*
* First spectrum is a detector. Remaining two spectra are monitors
* Detector contains histogram of tof with a Gaussian profile
*/
MatrixWorkspace_sptr CreateHistogramWorkspace()
{
const int numHists(3);
......@@ -138,14 +247,16 @@ MatrixWorkspace_sptr CreateHistogramWorkspace()
testWS->getAxis(0)->unit()=Mantid::Kernel::UnitFactory::Instance().create("TOF");
MantidVecPtr xdata;
xdata.access().resize(numBins+1);
const double peakHeight(1000), peakCentre(7000.), sigmaSq(1000*1000.);
const double peakHeight(1000.), peakCentre(7000.), sigmaSq(1000*1000.);
// tof ranges from 0 to 16000 (units assumed micro-seconds
const double rescaling_factor(4.0);
for(int ibin=0; ibin<numBins; ++ibin)
{
const double xValue=4*ibin;
const double xValue = rescaling_factor*ibin;
testWS->dataY(0)[ibin]=peakHeight*exp(-0.5*pow(xValue-peakCentre, 2.)/sigmaSq);
xdata.access()[ibin] = xValue;
}
xdata.access()[numBins] = 4*numBins;
xdata.access()[numBins] = rescaling_factor*numBins;
for( int ihist=0; ihist<numHists; ihist++)
testWS->setX(ihist, xdata);
return testWS;
......@@ -157,6 +268,7 @@ EventWorkspace_sptr CreateEventWorkspace()
const bool clearEvents(true);
EventWorkspace_sptr testWS=WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(numBanks,numPixels,clearEvents);
testWS->getAxis(0)->unit()=Mantid::Kernel::UnitFactory::Instance().create("TOF");
const double rescaling_factor(4.0);
const size_t numHists=testWS->getNumberHistograms();
for (size_t ihist=0; ihist<numHists; ++ihist)
{
......@@ -165,7 +277,7 @@ EventWorkspace_sptr CreateEventWorkspace()
xdata.access().resize(numBins+1);
for(int ibin=0; ibin<=numBins; ++ibin)
{
double tof=4*ibin;
double tof = rescaling_factor*ibin;
TofEvent tofevent(tof);
xdata.access()[ibin]=tof;
evlist.addEventQuickly(tofevent); // insert event
......@@ -175,16 +287,21 @@ EventWorkspace_sptr CreateEventWorkspace()
return testWS;
}
void AddToInstrument(MatrixWorkspace_sptr &testWS, const bool &add_deltaE_mode=false, const bool &add_t0_formula=false)
void AddToIndirectInstrument(MatrixWorkspace_sptr &testWS, const bool &add_t0_formula=false, const bool &add_Efixed=false)
{
const double evalue(2.082); // energy corresponding to the first order Bragg peak in the analyzers
if(add_deltaE_mode)
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),"deltaE-mode", "indirect");
for(size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
testWS->instrumentParameters().addDouble(testWS->getDetector(ihist)->getComponentID(),"Efixed",evalue);
if(add_t0_formula)
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),"t0_formula","50.-(50./52500)*incidentEnergy");
}
if(add_Efixed)
{
const double evalue(2.082); // energy corresponding to the first order Bragg peak in the analyzers
for(size_t ihist=0; ihist<testWS->getNumberHistograms(); ++ihist)
{
testWS->instrumentParameters().addDouble(testWS->getDetector(ihist)->getComponentID(),"Efixed",evalue);
}
}
} // end of void AddToInstrumen
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment