Skip to content

Commit

Permalink
Refs #6670 Finalize debugging of test file
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmborr
jmborr committed May 15, 2013
1 parent d1ac766 commit 36c8bb4
Show file tree
Hide file tree
Showing 2 changed files with 87 additions and 101 deletions.
39 changes: 21 additions & 18 deletions Code/Mantid/Framework/Algorithms/src/ModeratorTzero2.cpp
View file
Expand Up @@ -84,12 +84,12 @@ void ModeratorTzero2::exec()

//deltaE-mode (should be "indirect")
std::vector<std::string> Emode=m_instrument->getStringParameter("deltaE-mode");
if(!Emode.size()) throw Exception::InstrumentDefinitionError("Unable to retrieve instrument geometry (direct or indirect) parameter", inputWS->getTitle());
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.");

// extract formula from instrument parameters
std::vector<std::string> t0_formula=m_instrument->getStringParameter("t0_formula");
if(!t0_formula.size()) throw Exception::InstrumentDefinitionError("Unable to retrieve t0_formula among instrument parameters");
if(t0_formula.empty()) throw Exception::InstrumentDefinitionError("Unable to retrieve t0_formula among instrument parameters");
m_formula=t0_formula[0];

//Run execEvent if eventWorkSpace
Expand Down Expand Up @@ -125,14 +125,19 @@ void ModeratorTzero2::exec()
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++)
{
tof=inbins[ibin]; // current time-of-flight
if(tof<0.0) tof=0.0; // impossible substraction to negative time-of-flight values
outbins[ibin] = inbins[ibin]-CalculateT0(tof, L1, t2, E1, parser);
if(tof<m_t1min+t2)
tof-=min_t0_next;
else
tof-=CalculateT0(tof, L1, t2, E1, parser);
outbins[ibin] = tof;
}
}
else
Expand Down Expand Up @@ -216,35 +221,33 @@ void ModeratorTzero2::execEvent()
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1); // associate variable 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 parameter
// fix the histogram bins
MantidVec & x = evlist.dataX();
for (MantidVec::iterator iter = x.begin(); iter != x.end(); ++iter)
{
tof = *iter;
if(tof<0.0) tof=0.0; // impossible substraction to negative time-of-flights
*iter-=CalculateT0(tof, L1, t2, E1, parser);
if(tof<m_t1min+t2)
tof-=min_t0_next;
else
tof-=CalculateT0(tof, L1, t2, E1, parser);
*iter=tof;
}

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 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)
{
tofs[itof]=tof-min_t0_next;
}
tof-=min_t0_next;
else
{
tofs[itof]=tof-CalculateT0(tof, L1, t2, E1, parser);
}
tof-=CalculateT0(tof, L1, t2, E1, parser);
tofs[itof]=tof;
}
evlist.setTofs(tofs);
//evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
//evlist.sortTof();
evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
}
}
prog.report();
Expand Down
149 changes: 66 additions & 83 deletions Code/Mantid/Framework/Algorithms/test/ModeratorTzero2Test.h
View file
Expand Up @@ -3,6 +3,7 @@

#include <cxxtest/TestSuite.h>
#include "MantidAlgorithms/ModeratorTzero2.h"
#include "MantidAlgorithms/Rebin.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidDataObjects/EventList.h"
Expand All @@ -13,7 +14,6 @@
#include "MantidGeometry/muParser_Silent.h"
#include <cmath>


using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Algorithms;
Expand All @@ -36,7 +36,6 @@ void testExecThrowsDeltaEmode()
{
const std::string inputName("testWS");
Workspace2D_sptr testWS = createTestHistoWorkspace();
AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service
ModeratorTzero2 alg;
initializeAlgorithm(alg,inputName,inputName);
TS_ASSERT_THROWS(alg.execute(), Mantid::Kernel::Exception::InstrumentDefinitionError);
Expand All @@ -48,10 +47,7 @@ void testExecThrowsTzeroFormula()
{
Workspace2D_sptr testWS = createTestHistoWorkspace();
const std::string inputName("testWS");
AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service
Mantid::Geometry::Instrument_const_sptr instrument = testWS->getInstrument(); // pointer to the instrument
testWS->instrumentParameters().addString(instrument->getChild(0).get(),"deltaE-mode", "indirect");
std::vector<std::string> Emode=instrument->getStringParameter("deltaE-mode");
testWS->instrumentParameters().addString(testWS->getInstrument()->getComponentID(),"deltaE-mode", "indirect");
ModeratorTzero2 alg;
initializeAlgorithm(alg,inputName,inputName);
TS_ASSERT_THROWS(alg.execute(), Mantid::Kernel::Exception::InstrumentDefinitionError);
Expand All @@ -62,86 +58,59 @@ void testExecThrowsTzeroFormula()
void testExecHistogramWorkspace()
{
Workspace2D_sptr testWS = createTestHistoWorkspace();
const std::string inputName("testWS");
AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service
const std::string formula("0.1*incidentEnergy^(-0.25)");
const std::string formula("10*incidentEnergy^(-0.05)"); // formula for delay emission time from moderator versus incident energy
addFormula(formula, boost::dynamic_pointer_cast<MatrixWorkspace>(testWS)); // insert t0-formula as an instrument parameter
ModeratorTzero2 alg;
initializeAlgorithm(alg,inputName,inputName);
initializeAlgorithm(alg,"testWS","testWS");
TS_ASSERT_THROWS_NOTHING(alg.execute());
const int numHists(2);
const int numBins(2000);
double L1,L2,E2,v2,t2,tof,t1,v1,E1;
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
parser.SetExpr(formula);
Mantid::Geometry::IObjComponent_const_sptr sample = testWS->getInstrument()->getSample();
L1=testWS->getInstrument()->getSource()->getDistance(*sample);
for(int iHist=0; iHist<numHists; iHist++)
double tofs[4][11]={ // check a few values
{294.191, 1294.19, 2294.19, 3294.19, 4294.19, 5294.19, 6294.19, 7294.19, 8293.81, 9292.95, 10292.56},
{294.191, 1294.19, 2294.19, 3294.19, 4294.19, 5294.19, 6294.19, 7294.19, 8293.81, 9292.96, 10292.56},
{-6.16667, 193.833, 393.391, 593.117, 792.916, 992.757, 1192.62, 1392.51, 1592.41, 1792.32, 1992.24},
{-5.87034, 194.13, 393.708, 593.448, 793.257, 993.105, 1192.98, 1392.87, 1592.77, 1792.69, 1992.61}
};
for(unsigned int iHist=0; iHist<testWS->getNumberHistograms(); iHist++)
{
Mantid::MantidVec &xdata=testWS->dataX(iHist);
Mantid::Geometry::IDetector_const_sptr det=testWS->getDetector(iHist);
L2=det->getDistance(*sample);
std::vector<double> wsProp=det->getNumberParameter("Efixed");
E2=wsProp.at(0); //[E2]=meV
v2=m_convFact*sqrt(E2);
t2=v2/L2;
for(int iBin=0; iBin<=numBins; ++iBin)
{
tof = 5.0 + 5.5*iBin; // original time of flight
t1=tof-t2;
v1=L1/t1;
E1=m_convFact*v1*v1; // Energy in meV if v1 in meter/microsecond
TS_ASSERT_DELTA(xdata[iBin], tof-parser.Eval(), 1e-08); // evaluate the formula on the X-axis
}
for(unsigned int iBin=0; iBin<xdata.size(); iBin++)
if(!(iBin%200))
TS_ASSERT_DELTA(xdata[iBin], tofs[iHist][iBin/200], 0.01); // 0.01micro-second is more than enough precision
}
AnalysisDataService::Instance().remove(inputName); // unregister the workspace from the analysis data service
AnalysisDataService::Instance().remove("testWS"); // unregister the workspace from the analysis data service
}

/// Test transformation on an event workspace
void testExecEventWorkspace()
{
Mantid::DataObjects::EventWorkspace_sptr testWS=createTestEventWorkspace();
const std::string inputName("testWS");
AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service
const std::string formula("0.1*incidentEnergy^(-0.25)");
const std::string formula("10*incidentEnergy^(-0.05)"); // formula for delay emission time from moderator versus incident energy
addFormula(formula, boost::dynamic_pointer_cast<MatrixWorkspace>(testWS)); // insert t0-formula as an instrument parameter
ModeratorTzero2 alg;
const std::string outputName("outWS");
initializeAlgorithm(alg,inputName,outputName);
initializeAlgorithm(alg,"testWS","outWS");
TS_ASSERT_THROWS_NOTHING(alg.execute());
MatrixWorkspace_sptr outws = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>(outputName);
// check a few values
double tofH[17]={-6.22573, 993.774, 1993.77, 2993.77, 3993.77, 4993.77, 5993.77, 6993.77, 7993.77, 8992.63, 9992.13, 10991.8170, 11991.5835, 12991.4, 13991.2408, 14991.1, 15990.9881};
double tofL[10]={143.774, 593.774, 2993.77, 6993.77, 7993.77, 8992.63, 9992.13, 10991.8, 12991.4, 14991.1};
MatrixWorkspace_sptr outws = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("outWS");
Mantid::DataObjects::EventWorkspace_sptr outWS(boost::dynamic_pointer_cast<Mantid::DataObjects::EventWorkspace>(outws));
Mantid::Geometry::IObjComponent_const_sptr sample = testWS->getInstrument()->getSample();
double L1,L2,E2,v2,t2,t1,v1,E1;
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1); // associate E1 to this parser
parser.SetExpr(formula);
L1=testWS->getInstrument()->getSource()->getDistance(*sample);
const size_t numHists=testWS->getNumberHistograms();
for(size_t iHist=0; iHist<numHists; ++iHist)

for(size_t iHist=0; iHist<testWS->getNumberHistograms(); ++iHist)
{
Mantid::Geometry::IDetector_const_sptr det=testWS->getDetector(iHist);
L2=det->getDistance(*sample);
std::vector< double > wsProp=det->getNumberParameter("Efixed");
E2=wsProp.at(0); //[E2]=meV
v2=m_convFact*sqrt(E2);
t2=v2/L2;
Mantid::DataObjects::EventList &inevlist=testWS->getEventList(iHist);
Mantid::MantidVec intofs=inevlist.getTofs();
Mantid::DataObjects::EventList &outevlist=outWS->getEventList(iHist);
Mantid::MantidVec outtofs=outevlist.getTofs();
for(unsigned int itof=0; itof<intofs.size(); itof++)
Mantid::DataObjects::EventList &evlist=outWS->getEventList(iHist);
Mantid::MantidVec & x = evlist.dataX();
for(unsigned int iBin=0; iBin<x.size(); iBin++)
{
TS_ASSERT_DELTA(x[iBin], tofH[iBin], 0.01); // 0.01 micro-seconds accuracy for the bin boundaries
}
Mantid::MantidVec tofs=evlist.getTofs();
for(unsigned int itof=0; itof<tofs.size(); itof++)
{
t1=intofs[itof]-t2;
v1=L1/t1;
E1=m_convFact*v1*v1; // Energy in meV if v1 in meter/microsecond
TS_ASSERT_DELTA(outtofs[itof], intofs[itof]-parser.Eval(), 1e-08); //evaluate the formula on the X-axis
TS_ASSERT_DELTA(tofs[itof], tofL[itof], 0.2); //0.1 micro-seconds is the accuracy limit of this algorithm, thus 0.2 is a safe accuracy for testing
}

}
AnalysisDataService::Instance().remove(inputName); // unregister the workspace from the analysis data service
AnalysisDataService::Instance().remove(outputName); // unregister the workspace from the analysis data service
AnalysisDataService::Instance().remove("testWS"); // unregister the workspace from the analysis data service
AnalysisDataService::Instance().remove("outWS"); // unregister the workspace from the analysis data service
}

private:
Expand Down Expand Up @@ -169,45 +138,59 @@ Workspace2D_sptr createTestHistoWorkspace()
}
testWS->setX(0, xdata);
testWS->setX(1, xdata);
const std::string inputName("testWS");
AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service
return testWS;
}

/// Create an event workspace with time-of-flights values normally distributed
/// Create an event workspace containing two detectors and two monitors
Mantid::DataObjects::EventWorkspace_sptr createTestEventWorkspace()
{
const int numBanks(1); // one bank per instrument
const int numPixels(2); // Two pixels per bank
const int numEvents(200); // number of events per pixel
const int numBanks(1); // one squared bank per instrument
const int numPixels(2); // 2x2=4 pixels per bank
const int numEvents(10); // number of events per pixel
double tofs[10]={150.0, 600.0, 3000.0, 7000.0, 8000.0, 9000.0, 10000.0, 11000.0, 13000.0, 15000.0}; // use this list of TOF at each detector pixel
Mantid::DataObjects::EventWorkspace_sptr testWS=WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(numBanks, numPixels);
testWS->getAxis(0)->unit()=Mantid::Kernel::UnitFactory::Instance().create("TOF");
boost::mt19937 rng;
boost::normal_distribution<> nd(6493.0, 724.0); //normal distribution with mean=6493 and standard distribution=724
boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > var_nor(rng, nd); // random number generator
const size_t numHists=testWS->getNumberHistograms();
ModeratorTzero2 alg;
const double t1min(alg.gett1min());
double tof;
for(size_t iHist=0; iHist<numHists; ++iHist)
{
Mantid::DataObjects::EventList &evlist=testWS->getEventList(iHist);
for(int iEvent=0; iEvent<numEvents; ++iEvent)
{
tof=var_nor();
while(tof<t1min) tof=var_nor(); // fast neutrons are treated specially
Mantid::DataObjects::TofEvent event(tof);
Mantid::DataObjects::TofEvent event(tofs[iEvent]);
evlist.addEventQuickly(event);
}
}
const std::string inputName("testWS");
AnalysisDataService::Instance().add(inputName, testWS); // register the workspace in the analysis data service
Rebin alg; // Rebin the companion histogram
TS_ASSERT_THROWS_NOTHING(alg.initialize());
TS_ASSERT( alg.isInitialized() );
alg.setProperty("InputWorkspace",testWS);
alg.setPropertyValue("OutputWorkspace","testWS");
std::vector<double> params {0.,1000.,16000.};
alg.setProperty("Params",params);
alg.setProperty("PreserveEvents",true);
TS_ASSERT_THROWS_NOTHING(alg.execute());
return testWS;
}

/// Add formula for the emission time from moderator to instrument parameters
/* Add formula for the emission time from moderator to instrument parameters
* Also define instrument as of type "indirect"
* Also define "Efixed" for each detector (but nor for the monitors)
*/
void addFormula(const std::string &formula, MatrixWorkspace_sptr workspace)
{
//deltaE-mode (should be "indirect")
Mantid::Geometry::Instrument_const_sptr instrument = workspace->getInstrument(); // pointer to the instrument
workspace->instrumentParameters().addString(instrument->getChild(0).get(),"deltaE-mode", "indirect");
workspace->instrumentParameters().addString(instrument->getChild(0).get(),"t0_formula",formula);
const double efixed=2.082; // energy for analyzers, in meV
workspace->instrumentParameters().addString(workspace->getInstrument()->getComponentID(),"deltaE-mode", "indirect");
workspace->instrumentParameters().addString(workspace->getInstrument()->getComponentID(),"t0_formula", formula);
for (size_t i=0; i<workspace->getNumberHistograms(); ++i)
{
Mantid::Geometry::IDetector_const_sptr detector=workspace->getDetector(i);
if( !detector->isMonitor() )
workspace->instrumentParameters().addDouble(detector->getComponentID(),"Efixed", efixed);
}
}

/// Avoids writing these over and over
Expand Down

0 comments on commit 36c8bb4

Please sign in to comment.