Skip to content

Commit

Permalink
refs #5016. Extract and test q calculations
Browse files Browse the repository at this point in the history
  • Loading branch information
OwenArnold committed May 31, 2012
1 parent 193208f commit bd49912
Show file tree
Hide file tree
Showing 4 changed files with 165 additions and 27 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ namespace MDEvents
class DLLExport ReflectometryMDTransform
{
public:
virtual Mantid::API::IMDEventWorkspace_sptr execute(Mantid::API::IEventWorkspace_const_sptr eventWs) const = 0;
virtual Mantid::API::IMDEventWorkspace_sptr execute(Mantid::API::MatrixWorkspace_const_sptr inputWs) const = 0;
};
}
}
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -7,8 +7,94 @@

namespace Mantid
{
namespace MDEvents
{
namespace MDEvents
{
/**
Base class for reflectometry Q transformations
*/
class CalculateReflectometryQBase
{
protected:
const double to_radians_factor;
const double two_pi;
CalculateReflectometryQBase() : to_radians_factor(3.14159265/180), two_pi(6.28318531)
{
}
protected:
~CalculateReflectometryQBase(){};
};

/**
@class Converts from inputs of wavelength, incident theta and final theta to Qx for reflectometry experiments
*/
class CalculateReflectometryQx : public CalculateReflectometryQBase
{
private:
double m_cos_theta_i;
double m_dirQx;
public:
/**
Constructor
@param thetaIncident: incident theta value in degrees
*/
CalculateReflectometryQx(const double& thetaIncident): m_cos_theta_i(cos(thetaIncident*to_radians_factor))
{
}
/**
Setter for the final theta value require for the calculation. Internally pre-calculates and caches to cos theta for speed.
@param thetaFinal: final theta value in degrees
*/
void setThetaFinal(const double& thetaFinal)
{
const double c_cos_theta_f = cos(thetaFinal*to_radians_factor);
m_dirQx = (c_cos_theta_f - m_cos_theta_i);
}
/**
Executes the calculation to determine Qz
@param wavelength : wavelenght in Anstroms
*/
double execute(const double& wavelength) const
{
double wavenumber = two_pi/wavelength;
return wavenumber * m_dirQx;
}
};

/**
@class Converts from inputs of wavelength, incident theta and final theta to Qz for reflectometry experiments
*/
class CalculateReflectometryQz : public CalculateReflectometryQBase
{
private:
double m_sin_theta_i;
double m_dirQz;
public:
/**
Constructor
@param thetaIncident: incident theta value in degrees
*/
CalculateReflectometryQz(const double& thetaIncident): m_sin_theta_i(sin(thetaIncident*to_radians_factor))
{
}
/**
Setter for the final theta value require for the calculation. Internally pre-calculates and caches to sine theta for speed.
@param thetaFinal: final theta value in degrees
*/
void setThetaFinal(const double& thetaFinal)
{
const double c_sin_theta_f = sin(thetaFinal*to_radians_factor);
m_dirQz = (c_sin_theta_f + m_sin_theta_i);
}
/**
Executes the calculation to determine Qz
@param wavelength : wavelenght in Anstroms
*/
double execute(const double& wavelength) const
{
double wavenumber = two_pi/wavelength;
return wavenumber * m_dirQz;
}
};

/** ReflectometryTranformQxQz : Type of ReflectometyTransform. Used to convert from an input event workspace to a 2D MDEvent workspace with dimensions of QxQy.
Transformation is specific for reflectometry purposes.
Expand Down Expand Up @@ -42,15 +128,18 @@ namespace MDEvents
const double m_qxMax;
const double m_qzMin;
const double m_qzMax;
const double m_incidentTheta;
/// Object performing raw caclcation to determine Qx
mutable CalculateReflectometryQx m_QxCalculation;
/// Object performing raw calculation to determine Qx
mutable CalculateReflectometryQz m_QzCalculation;
public:

/// Constructor
ReflectometryTransformQxQz(double qxMin, double qxMax, double qzMin, double qzMax, double incidentTheta);
/// Destructor
~ReflectometryTransformQxQz();
/// Execute transformation
virtual Mantid::API::IMDEventWorkspace_sptr execute(Mantid::API::IEventWorkspace_const_sptr eventWs) const;
virtual Mantid::API::IMDEventWorkspace_sptr execute(Mantid::API::MatrixWorkspace_const_sptr inputWs) const;

private:

Expand Down
36 changes: 14 additions & 22 deletions Code/Mantid/Framework/MDEvents/src/ReflectometryTransformQxQz.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,7 @@ namespace Mantid
@param incidentTheta: Predetermined incident theta value
*/
ReflectometryTransformQxQz::ReflectometryTransformQxQz(double qxMin, double qxMax, double qzMin, double qzMax, double incidentTheta):
m_qxMin(qxMin), m_qxMax(qxMax), m_qzMin(qzMin), m_qzMax(qzMax), m_incidentTheta(incidentTheta)
m_qxMin(qxMin), m_qxMax(qxMax), m_qzMin(qzMin), m_qzMax(qzMax), m_QxCalculation(incidentTheta), m_QzCalculation(incidentTheta)
{
if(qxMin >= qxMax)
{
Expand All @@ -47,9 +47,9 @@ namespace Mantid
/*
Execute the transformtion. Generates an output IMDEventWorkspace.
@return the constructed IMDEventWorkspace following the transformation.
@param ws: Input EventWorkspace const shared pointer
@param ws: Input MatrixWorkspace const shared pointer
*/
IMDEventWorkspace_sptr ReflectometryTransformQxQz::execute(IEventWorkspace_const_sptr eventWs) const
IMDEventWorkspace_sptr ReflectometryTransformQxQz::execute(MatrixWorkspace_const_sptr inputWs) const
{
const size_t nbinsx = 10;
const size_t nbinsz = 10;
Expand All @@ -72,30 +72,22 @@ namespace Mantid
// Start with a MDGridBox.
ws->splitBox();

auto spectraAxis = eventWs->getAxis(1);
const double two_pi = 6.28318531;
const double to_radians_factor = 3.14159265/180;
const double c_cos_theta_i = cos(m_incidentTheta*to_radians_factor);
const double c_sin_theta_i = sin(m_incidentTheta*to_radians_factor);
for(size_t index = 0; index < eventWs->getNumberHistograms(); ++index)
auto spectraAxis = inputWs->getAxis(1);
for(size_t index = 0; index < inputWs->getNumberHistograms(); ++index)
{
auto counts = eventWs->readY(index);
auto wavelengths = eventWs->readX(index);
auto errors = eventWs->readE(index);
size_t nInputBins = eventWs->isHistogramData() ? wavelengths.size() -1 : wavelengths.size();
auto counts = inputWs->readY(index);
auto wavelengths = inputWs->readX(index);
auto errors = inputWs->readE(index);
size_t nInputBins = inputWs->isHistogramData() ? wavelengths.size() -1 : wavelengths.size();
const double theta_final = spectraAxis->getValue(index);
const double c_sin_theta_f = sin(theta_final*to_radians_factor);
const double c_cos_theta_f = cos(theta_final*to_radians_factor);
const double dirQx = (c_cos_theta_f - c_cos_theta_i);
const double dirQz = (c_sin_theta_f + c_sin_theta_i);
m_QxCalculation.setThetaFinal(theta_final);
m_QzCalculation.setThetaFinal(theta_final);
//Loop over all bins in spectra
for(size_t binIndex = 0; binIndex < nInputBins; ++binIndex)
{
double lambda = wavelengths[binIndex];
double wavenumber = two_pi/lambda;
double _qx = wavenumber * dirQx;
double _qz = wavenumber * dirQz;

const double& lambda = wavelengths[binIndex];
double _qx = m_QxCalculation.execute(lambda);
double _qz = m_QzCalculation.execute(lambda);
double centers[2] = {_qx, _qz};

ws->addEvent(MDLeanEvent<2>(float(counts[binIndex]), float(errors[binIndex]*errors[binIndex]), centers));
Expand Down
Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
#ifndef MANTID_MDEVENTS_REFLECTOMETRYTRANFORMQXQZTEST_H_
#define MANTID_MDEVENTS_REFLECTOMETRYTRANFORMQXQZTEST_H_

#define PI 3.14159

#include <cxxtest/TestSuite.h>
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
Expand All @@ -15,6 +17,8 @@ using namespace Mantid::API;

class ReflectometryTransformQxQzTest : public CxxTest::TestSuite
{
private:

public:
// This pair of boilerplate methods prevent the suite being created statically
// This means the constructor isn't called when running other tests
Expand Down Expand Up @@ -92,6 +96,59 @@ class ReflectometryTransformQxQzTest : public CxxTest::TestSuite
TS_ASSERT_THROWS_NOTHING(ReflectometryTransformQxQz(qxMin, qxMax, qzMin, qzMax, incidentTheta));
}

//---- Tests for Qx Calculator ---- //

void test_calculate_Qx()
{
//Set up calculation so that it collapses down to 2*PI/wavelength by setting initial theta to PI/2 and final theta to zero
CalculateReflectometryQx calculator(90);
double qx;
const double wavelength = 0.1;
TS_ASSERT_THROWS_NOTHING(calculator.setThetaFinal(0));
TS_ASSERT_THROWS_NOTHING(qx = calculator.execute(wavelength));
TS_ASSERT_DELTA(2*PI/wavelength, qx, 0.0001);
}

void test_recalculate_Qx()
{
CalculateReflectometryQx calculator(0);
calculator.setThetaFinal(0);
const double wavelength = 0.1;
TS_ASSERT_DELTA(0, calculator.execute(wavelength), 0.0001);

//Now reset the final theta and should be able to re-execute
calculator.setThetaFinal(90);
TS_ASSERT_DELTA(-2*PI/wavelength, calculator.execute(wavelength), 0.0001);
}

//---- End Tests for Qx Calculator ---- //

//---- Tests for Qz Calculator ---- //

void test_calculate_Qz()
{
//Set up calculation so that it collapses down to 2*PI/wavelength
CalculateReflectometryQz calculator(0);
double qx;
const double wavelength = 0.1;
TS_ASSERT_THROWS_NOTHING(calculator.setThetaFinal(90));
TS_ASSERT_THROWS_NOTHING(qx = calculator.execute(wavelength));
TS_ASSERT_DELTA(2*PI/wavelength, qx, 0.0001);
}

void test_recalculate_Qz()
{
CalculateReflectometryQz calculator(90);
calculator.setThetaFinal(90);
const double wavelength = 0.1;
TS_ASSERT_DELTA(2*(2*PI/wavelength), calculator.execute(wavelength), 0.001);

//Now reset the final theta and should be able to re-execute
calculator.setThetaFinal(0);
TS_ASSERT_DELTA(2*PI/wavelength, calculator.execute(wavelength), 0.001);
}

//---- End Tests for Qz Calculator ---- //
};


Expand Down

0 comments on commit bd49912

Please sign in to comment.