Skip to content

Commit

Permalink
Refs #6670 Fixed problems with unphysical events
Browse files Browse the repository at this point in the history
  • Loading branch information
@jmborr
jmborr committed May 2, 2013
1 parent 9b83954 commit 639792a
Show file tree
Hide file tree
Showing 3 changed files with 105 additions and 62 deletions.
7 changes: 5 additions & 2 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ModeratorTzero2.h
View file
Expand Up @@ -5,6 +5,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/Algorithm.h"
#include "MantidKernel/PhysicalConstants.h"

namespace Mantid
{
Expand Down Expand Up @@ -61,7 +62,7 @@ class DLLExport ModeratorTzero2 : public Mantid::API::Algorithm
{
public:
/// (Empty) Constructor
ModeratorTzero2() : Mantid::API::Algorithm(), m_niter(10) {}
ModeratorTzero2() : Mantid::API::Algorithm(), m_convfactor(0.5e+12*Mantid::PhysicalConstants::NeutronMass/Mantid::PhysicalConstants::meV), m_niter(1) {}
/// Virtual destructor
virtual ~ModeratorTzero2() {}
/// Algorithm's name
Expand Down Expand Up @@ -89,8 +90,10 @@ class DLLExport ModeratorTzero2 : public Mantid::API::Algorithm
double CalculateE1(MantidVec &tofs, unsigned int j, double t0, double t2, const double &L1);
/// Calculate the emission time from the moderator
double EvaluateFormula(const std::string &formula, const double &energy);

const double m_convfactor;
/// Maximum number of iterations when calculating the emission time from the moderator
const size_t m_niter;
size_t m_niter;
/// tolerance for calculating E1, in picoseconds
double m_tolTOF;
/// string containing the heuristic regression for the moderator emission time versus neutron energy
Expand Down
151 changes: 92 additions & 59 deletions Code/Mantid/Framework/Algorithms/src/ModeratorTzero2.cpp
View file
Expand Up @@ -34,6 +34,7 @@
#include "MantidKernel/UnitFactory.h"
#include "MantidGeometry/muParser_Silent.h"
#include <boost/lexical_cast.hpp>
#include "MantidDataObjects/EventList.h"

namespace Mantid
{
Expand Down Expand Up @@ -66,37 +67,15 @@ void ModeratorTzero2::init()
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");
//declare the output workspace
declareProperty(new WorkspaceProperty<MatrixWorkspace>("OutputWorkspace","",Direction::Output), "The name of the output workspace");
declareProperty(new Kernel::PropertyWithValue<double>("tolTOF", 4.0, Kernel::Direction::Input),"Tolerance in the calculation of the emission time, in microseconds");
declareProperty(new Kernel::PropertyWithValue<double>("tolTOF", 1.0, Kernel::Direction::Input),"Tolerance in the calculation of the emission time, in microseconds");
declareProperty(new Kernel::PropertyWithValue<size_t>("niter", 1, Kernel::Direction::Input),"Number of iterations (default:1)");

} // end of void ModeratorTzero2::init()

double ModeratorTzero2::EvaluateFormula(const std::string &formula, const double &energy)
{
std::string fe(formula);
std::stringstream e;
e<<energy;
size_t found;
while ( fe.find("incidentEnergy") != std::string::npos )
{
found = fe.find("incidentEnergy");
fe.replace(found, 14, e.str());
}
try
{
mu::Parser p;
p.SetExpr(fe);
return p.Eval();
}
catch (mu::Parser::exception_type &e)
{
throw Kernel::Exception::InstrumentDefinitionError("Equation attribute for t0 formula. Muparser error message is: " + e.GetMsg());
}
return 0;
}

void ModeratorTzero2::exec()
{
m_tolTOF = getProperty("tolTOF"); //time unit increment, in picoseconds;
m_niter=getProperty("niter"); // number of iterations
const MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
m_instrument = inputWS->getInstrument(); // pointer to the instrument

Expand All @@ -120,13 +99,12 @@ void ModeratorTzero2::exec()
// extract formula from instrument parameters
const ParameterMap& pmap = inputWS->constInstrumentParameters();
Instrument_const_sptr instrument=inputWS->getInstrument();
Parameter_sptr par = pmap.getRecursive(instrument->getChild(0).get(),"t0_formula");
Parameter_sptr par=pmap.getRecursive(instrument->getChild(0).get(),"t0_formula");
if (!par)
{
throw std::invalid_argument("t0_formula not found in parameters");
}
m_formula=par->asString();

//Run execEvent if eventWorkSpace
EventWorkspace_const_sptr eventWS = boost::dynamic_pointer_cast<const EventWorkspace>(inputWS);
if (eventWS != NULL)
Expand All @@ -144,36 +122,43 @@ void ModeratorTzero2::exec()
}

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)
//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)
{
PARALLEL_START_INTERUPT_REGION

//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
if(t2 >= 0) //t2 < 0 when no detector info is available
{
double tof, t0_curr, t0_next, t1, v1, E1; // local variables
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1);
parser.SetExpr(m_formula);
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 t0;
double t0prime(0.0); // initial emission time set to zero
tof=inbins[ibin]; // current time-of-flight
if(tof<0.0) tof=0.0; // impossible substraction to negative time-of-flights
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(t0prime-t0)>m_tolTOF && iiter<m_niter)
while (std::fabs(t0_curr-t0_next)>=m_tolTOF && iiter<m_niter)
{
t0=t0prime;
double E1=CalculateE1(inbins, ibin, t0, t2, L1);
t0prime=EvaluateFormula(m_formula, E1);
t0_curr=t0_next;
t1=tof-t0_curr-t2;
v1=L1/t1;
E1=m_convfactor*v1*v1; // Energy in meV if v1 in meter/microsecond
t0_next=parser.Eval();
iiter++;
}
outbins[ibin] = inbins[ibin] - t0;
outbins[ibin] = inbins[ibin] - t0_next;
}
}
else
Expand All @@ -183,10 +168,10 @@ void ModeratorTzero2::exec()
//Copy y and e data
outputWS->dataY(i) = inputWS->dataY(i);
outputWS->dataE(i) = inputWS->dataE(i);
prog.report();
PARALLEL_END_INTERUPT_REGION
//prog.report();
//PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
//PARALLEL_CHECK_INTERUPT_REGION

// Copy units
if (inputWS->getAxis(0)->unit().get())
Expand Down Expand Up @@ -218,11 +203,10 @@ void ModeratorTzero2::exec()
*/
double ModeratorTzero2::CalculateE1(MantidVec &tofs, unsigned int ibin, double t0, double t2, const double &L1)
{
static const double convfactor = 0.5e+12*PhysicalConstants::NeutronMass/PhysicalConstants::meV;
double E1(0.0);
double tof=tofs[ibin];
double v1 = L1/(tof-t0-t2);
E1=convfactor*v1*v1; // Energy in meV if v1 in meter/microsecond
E1=m_convfactor*v1*v1; // Energy in meV if v1 in meter/microsecond
return E1;
}

Expand Down Expand Up @@ -258,41 +242,90 @@ void ModeratorTzero2::execEvent()
IObjComponent_const_sptr sample = outputWS->getInstrument()->getSample();

// Loop over the spectra
Progress prog(this,0.0,1.0,numHists); //report progress of algorithm
PARALLEL_FOR1(outputWS)
//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)
{
//PARALLEL_START_INTERUPT_REGION
size_t wsIndex = static_cast<size_t>(i);
PARALLEL_START_INTERUPT_REGION
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
{
MantidVec tofs=evlist.getTofs();
for(unsigned int ibin=0; ibin<tofs.size(); ibin++) {
double t0;
double t0prime(0.0);
double tof, t0_curr, t0_next, t1, v1, E1; // local variables
mu::Parser parser;
parser.DefineVar("incidentEnergy", &E1);
parser.SetExpr(m_formula);

// fix the histogram parameter
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
t0_curr=m_tolTOF;
t0_next=0.0;
size_t iiter(0);
while(std::fabs(t0_curr-t0_next)>=m_tolTOF && iiter<m_niter)
{
t0_curr=t0_next;
t1=tof-t0_curr-t2;
v1=L1/t1;
E1=m_convfactor*v1*v1; // Energy in meV if v1 in meter/microsecond
t0_next=parser.Eval();
iiter++;
}
*iter-=t0_next;
}

t1=200.0; // fast neutrons are shifted by same amount
v1=L1/t1;
E1=m_convfactor*v1*v1;
double min_t0_next=parser.Eval();

MantidVec tofs=evlist.getTofs();
for(unsigned int itof=0; itof<tofs.size(); itof++) {
tof=tofs[itof];
t0_curr=m_tolTOF;
t0_next=0.0;
size_t iiter(0);
while(std::fabs(t0prime-t0)>m_tolTOF && iiter<m_niter)
while(std::fabs(t0_curr-t0_next)>=m_tolTOF && iiter<m_niter)
{
t0=t0prime;
double E1=CalculateE1(tofs, ibin, t0, t2, L1);
t0prime=EvaluateFormula(m_formula,E1);
t0_curr=t0_next;
t1=tof-t0_curr-t2;
if(t1<=200.0)
{
t0_next=min_t0_next;
break;
}
v1=L1/t1;
E1=m_convfactor*v1*v1; // Energy in meV if v1 in meter/microsecond
t0_next=parser.Eval();
iiter++;
}
tofs[ibin]-=t0; // remove the emission time
/*
if(tof-t0_next>=814.0 && tof-t0_next<=815.0)
{
//#tof t1 t2 t0_next tof-t0_next v1 E1 L1
std::cout<<tof<<" "<<t1<<" "<<t2<<" "<<t0_next<<" "<<tof-t0_next<<" "<<v1<<" "<<E1<<" "<<L1<<std::endl;
}
*/
std::cout<<tof<<" "<<t1<<" "<<t2<<" "<<t0_next<<" "<<tof-t0_next<<" "<<v1<<" "<<E1<<" "<<L1<<std::endl;
tofs[itof]-=t0_next; // remove the emission time
}
evlist.setTofs(tofs);
evlist.setSortOrder(Mantid::DataObjects::EventSortType::UNSORTED);
evlist.sortTof();
}
}
prog.report();
PARALLEL_END_INTERUPT_REGION
//prog.report();
//PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
outputWS->clearMRU(); // Clears the Most Recent Used lists */
//PARALLEL_CHECK_INTERUPT_REGION
outputWS->clearMRU(); // Clears the Most Recent Used lists */
} // end of void ModeratorTzero2::execEvent()

/// calculate distance from source to sample or detector
Expand Down
9 changes: 8 additions & 1 deletion Code/Mantid/instrument/VISION_Parameters.xml
View file
Expand Up @@ -60,7 +60,14 @@
<parameter name="Workflow.Masking" type="string">
<value val="IdentifyNoisyDetectors" />
</parameter>


<!-- Emission time from moderator, syntax below follows muParser, less-than sign escaped with the amperstand for correct xml parsing, incidentEnergy in meV, results in micro-seconds -->
<parameter name="t0_formula" type="string">
<!-- value val="(incidentEnergy &lt; 34.7332) ? 37.011296*incidentEnergy^(-0.052874) : (incidentEnergy &lt; 88.7556) ? 124.267307*incidentEnergy^(-0.394282) : (incidentEnergy &lt; 252.471) ? 963.775145*incidentEnergy^(-0.850919) : (incidentEnergy &lt; 420.145) ? 33.225834*incidentEnergy^(-0.242105) : (incidentEnergy &lt; 100000.0) ? 120.569231*incidentEnergy^(-0.455477) : 0.0" / -->
<value val="370*incidentEnergy^(-0.1)" />
<!-- value val="40.0"/-->
</parameter>

</component-link>

</parameter-file>

0 comments on commit 639792a

Please sign in to comment.