Skip to content

Commit

Permalink
refs #5016. Characterise current behaviour.
Browse files Browse the repository at this point in the history
Pause while we consider whether its possible or better to implement this algorithm in terms of others.
  • Loading branch information
OwenArnold committed Mar 23, 2012
1 parent f37cd1c commit 8aa1b08
Show file tree
Hide file tree
Showing 2 changed files with 5 additions and 128 deletions.
96 changes: 2 additions & 94 deletions Code/Mantid/Framework/MDEvents/src/ConvertToReflectometryQ.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -80,100 +80,8 @@ namespace MDEvents
*/
void ConvertToReflectometryQ::exec()
{

//TODO rebin input?

MatrixWorkspace_sptr in_ws = getProperty("InputWorkspace");

//TODO, how to calculate max, min + number of bins? Jon Taylor thinks that there is a way to calculate these upfront.
const double qzmin = -0.01;
const double qxmin = -0.01;
const double qxmax = 0.01;
const double qzmax = 0.01;
const size_t nbinsx = 10;
const size_t nbinsz = 10;

// Fixed dimensionality
MDHistoDimension_sptr qxDim = MDHistoDimension_sptr(new MDHistoDimension("Qx","qx","(Ang^-1)", qxmin, qxmax, nbinsx));
MDHistoDimension_sptr qzDim = MDHistoDimension_sptr(new MDHistoDimension("Qz","qz","(Ang^-1)", qzmin, qzmax, nbinsz));
MDHistoWorkspace_sptr ws = MDHistoWorkspace_sptr(new MDHistoWorkspace(qxDim, qzDim));

// check the input units
if (in_ws->getAxis(0)->unit()->unitID() != "TOF")
throw std::invalid_argument("Input event workspace's X axis must be in TOF units.");

// Calculate theta incident
Instrument_const_sptr instrument = in_ws->getInstrument();
V3D sample = instrument->getSample()->getPos();
V3D source = instrument->getSource()->getPos();
V3D beamDir = (sample - source);
V3D normBeamDir = beamDir / beamDir.norm();
V3D along = instrument->getReferenceFrame()->vecPointingAlongBeam();
V3D up = instrument->getReferenceFrame()->vecPointingUp();
double thetaIn = beamDir.angle(along);
double thetaFinal = 0;

const size_t nHistograms = in_ws->getNumberHistograms();
//Loop through all spectra
for(size_t index = 0; index < nHistograms; ++index)
{
V3D sink = in_ws->getDetector(index)->getPos();
V3D detectorDir = (sink - sample);
double distance = detectorDir.norm() + beamDir.norm();
V3D normalDetectorDir = detectorDir / detectorDir.norm();
thetaFinal = detectorDir.angle(along);

auto counts = in_ws->readY(index);
auto tofs = in_ws->readX(index);
auto errors = in_ws->readE(index);

size_t nTofs = in_ws->isHistogramData() ? (tofs.size() - 1) : tofs.size();

const double wavenumber_in_angstrom_times_tof_in_microsec =
(PhysicalConstants::NeutronMass * distance * 1e-10) / (1e-6 * PhysicalConstants::h_bar);

//Loop through all TOF
for(size_t tof = 0; tof < nTofs; ++tof)
{
//Calculate wave number
double t = tofs[tof];
if(in_ws->isHistogramData())
{
t = (t + tofs[tof+1])/2;
}
double wavenumber = wavenumber_in_angstrom_times_tof_in_microsec / t;

//ki - kf
V3D QDir = (normBeamDir - normalDetectorDir);
double _qx = QDir.X() * wavenumber;
double _qz = QDir.Z() * wavenumber;

/// If q-max and min are known a-prori, these boundrary case truncations are not required. See top of method for comment.
_qx = _qx < qxmin ? qxmin : _qx;
_qz = _qz < qzmin ? qzmin : _qz;
_qx = _qx > qxmax ? qxmax : _qx;
_qz = _qz > qzmax ? qzmax : _qz;

//Set up for linear transformation qi -> dimension index.
double mx = (nbinsx / (qxmax - qxmin));
double mz = (nbinsz / (qzmax - qzmin));
double cx = (nbinsx - mx * (qxmin + qxmax))/2;
double cz = (nbinsz - mz * (qzmin + qzmax))/2;

double posIndexX = mx*_qx + cx;
double posIndexZ = mz*_qz + cz;

size_t linearindex = (posIndexX * nbinsx) + posIndexZ;

double error = errors[tof];

//Do we want to accumulate signal values, i.e will there be any overlap in Q?
ws->setSignalAt(linearindex, ws->getSignalAt(linearindex) + counts[tof]);
ws->setErrorSquaredAt(linearindex, error*error);
}
}

setProperty("OutputWorkspace", ws);
//Impelementation details still to workout here. May be best/possible to implement this algorithm in terms of other algorithms.
throw std::runtime_error("Not Yet Implemented.");
}


Expand Down
37 changes: 3 additions & 34 deletions Code/Mantid/Framework/MDEvents/test/ConvertToReflectometryQTest.h
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,12 @@ class ConvertToReflectometryQTest : public CxxTest::TestSuite
in_ws->getAxis(0)->setUnit("TOF");

ConvertToReflectometryQ alg;
alg.setRethrows(true);
TS_ASSERT_THROWS_NOTHING( alg.initialize() )
TS_ASSERT( alg.isInitialized() )
TS_ASSERT_THROWS_NOTHING( alg.setProperty("InputWorkspace", in_ws) );
TS_ASSERT_THROWS_NOTHING( alg.setPropertyValue("OutputWorkspace", outWSName) );
TS_ASSERT_THROWS_NOTHING( alg.execute(); );
TS_ASSERT( alg.isExecuted() );
TSM_ASSERT_THROWS("Not Implemented Yet", alg.execute(), std::runtime_error );
}

public:
Expand All @@ -49,42 +49,11 @@ class ConvertToReflectometryQTest : public CxxTest::TestSuite
TS_ASSERT( alg.isInitialized() )
}

void test_check_output_metadata()
void test_check_shouldnt_be_working()
{
std::string outWSName("ConvertToReflectometryQTest_OutputWS");

doExecute(outWSName);

IMDHistoWorkspace_sptr ws;
TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<IMDHistoWorkspace>(outWSName) );
TS_ASSERT(ws);

//Quick metadata checks associated with the ouput workspace
TSM_ASSERT_EQUALS("Wrong number of bins", 100, ws->getNPoints());
TSM_ASSERT_EQUALS("Wrong number of dimensions", 2, ws->getNumDims());
TSM_ASSERT_EQUALS("Wrong dimension names", "Qx", ws->getDimension(0)->getName());
TSM_ASSERT_EQUALS("Wrong dimension names", "Qz", ws->getDimension(1)->getName());
TSM_ASSERT_EQUALS("Wrong units on qx", "(Ang^-1)", ws->getDimension(0)->getUnits());
TSM_ASSERT_EQUALS("Wrong units on qz", "(Ang^-1)", ws->getDimension(1)->getUnits());
TSM_ASSERT_EQUALS("Wrong number of bins in qx", 10, ws->getDimension(0)->getNBins());
TSM_ASSERT_EQUALS("Wrong number of bins in qy", 10, ws->getDimension(1)->getNBins());

AnalysisDataService::Instance().remove(outWSName);
}

void test_check_output_data()
{
std::string outWSName("ConvertToReflectometryQTest_OutputWS");

doExecute(outWSName);

IMDHistoWorkspace_sptr ws;
TS_ASSERT_THROWS_NOTHING( ws = AnalysisDataService::Instance().retrieveWS<IMDHistoWorkspace>(outWSName) );
TS_ASSERT(ws);

//Detailed checks associated with the output workspace.

AnalysisDataService::Instance().remove(outWSName);
}

};
Expand Down

0 comments on commit 8aa1b08

Please sign in to comment.