Skip to content

Commit

Permalink
refs #9184. Test point detector position correction.
Browse files Browse the repository at this point in the history
  • Loading branch information
@OwenArnold
OwenArnold committed Jun 2, 2014
1 parent 17c8685 commit d0d12dd
Show file tree
Hide file tree
Showing 3 changed files with 153 additions and 67 deletions.
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ReflectometryReductionOne.h
View file
Expand Up @@ -72,7 +72,7 @@ namespace Mantid

/// Correct detector positions.
void correctPosition(API::MatrixWorkspace_sptr toCorrect, const double& thetaInDeg,
Geometry::IComponent_const_sptr sample, Geometry::IComponent_const_sptr detector);
Geometry::IComponent_const_sptr sample, Geometry::IComponent_const_sptr detector, const bool isPointDetector);

/// Sum spectra.
Mantid::API::MatrixWorkspace_sptr sumSpectraOverRange(API::MatrixWorkspace_sptr inWS, const int startIndex, const int endIndex);
Expand Down
194 changes: 130 additions & 64 deletions Code/Mantid/Framework/Algorithms/src/ReflectometryReductionOne.cpp
View file
@@ -1,26 +1,26 @@
/*WIKI*
Reduces a single TOF reflectometry run into a mod Q vs I/I0 workspace. Performs transmission corrections. Handles both point detector and multidetector cases.
The algorithm can correct detector locations based on an input theta value.
Reduces a single TOF reflectometry run into a mod Q vs I/I0 workspace. Performs transmission corrections. Handles both point detector and multidetector cases.
The algorithm can correct detector locations based on an input theta value.
Historically the work performed by this algorithm was known as the Quick script.
Historically the work performed by this algorithm was known as the Quick script.
=== Analysis Modes ===
=== Analysis Modes ===
The default analysis mode is ''PointDetectorAnalysis''. Only this mode supports Transmission corrections (see below). For PointAnalysisMode
the analysis can be roughly reduced to IvsLam = DetectorWS / sum(I0) / TransmissionWS / sum(I0). The normalization by tranmission run(s) is optional.
Input workspaces are converted to ''Wavelength'' first via [[ConvertUnits]].
The default analysis mode is ''PointDetectorAnalysis''. Only this mode supports Transmission corrections (see below). For PointAnalysisMode
the analysis can be roughly reduced to IvsLam = DetectorWS / sum(I0) / TransmissionWS / sum(I0). The normalization by tranmission run(s) is optional.
Input workspaces are converted to ''Wavelength'' first via [[ConvertUnits]].
IvsQ is calculated via [[ConvertUnits]] into units of ''MomentumTransfer''.
Corrections may be applied prior to the transformation to ensure that the detectors are in the correct location according to the input Theta value.
Corrections are only enabled when a Theta input value has been provided.
IvsQ is calculated via [[ConvertUnits]] into units of ''MomentumTransfer''.
Corrections may be applied prior to the transformation to ensure that the detectors are in the correct location according to the input Theta value.
Corrections are only enabled when a Theta input value has been provided.
=== Transmission Runs ===
Transmission correction is a normalization step, which may be applied to ''PointDetectorAnalysis'' reduction.
=== Transmission Runs ===
Transmission correction is a normalization step, which may be applied to ''PointDetectorAnalysis'' reduction.
Transmission runs are expected to be in TOF. The spectra numbers in the Transmission run workspaces must be the same as those in the Input Run workspace.
If two Transmission runs are provided then the Stitching parameters associated with the transmission runs will also be required.
If a single Transmission run is provided, then no stitching
parameters will be needed.
Transmission runs are expected to be in TOF. The spectra numbers in the Transmission run workspaces must be the same as those in the Input Run workspace.
If two Transmission runs are provided then the Stitching parameters associated with the transmission runs will also be required.
If a single Transmission run is provided, then no stitching
parameters will be needed.
*WIKI*/

Expand Down Expand Up @@ -79,6 +79,25 @@ namespace Mantid
const std::string pointDetectorAnalysis = "PointDetectorAnalysis";
const std::string tofUnitId = "TOF";
const std::string wavelengthUnitId = "Wavelength";

/**
* Helper free function to get the ordered spectrum numbers from a workspace.
* @param ws
* @return
*/
std::vector<int> getSpectrumNumbers(MatrixWorkspace_sptr& ws)
{
auto specToWSIndexMap = ws->getSpectrumToWorkspaceIndexMap();
std::vector<int> keys(specToWSIndexMap.size());
size_t i = 0;
for (auto it = specToWSIndexMap.begin(); it != specToWSIndexMap.end(); ++it, ++i)
{
keys[i] = static_cast<int>(it->first);
++i;
}
std::sort(keys.begin(), keys.end()); // Sort the keys, as the order is not guaranteed in the map.
return keys;
}
}
/* End of ananomous namespace */

Expand Down Expand Up @@ -175,17 +194,16 @@ namespace Mantid
declareProperty(
new WorkspaceProperty<MatrixWorkspace>("FirstTransmissionRun", "", Direction::Input,
PropertyMode::Optional),
"First transmission run, or the low wavelength transmision run if SecondTransmissionRun is also provided.");
"First transmission run, or the low wavelength transmission run if SecondTransmissionRun is also provided.");
declareProperty(
new WorkspaceProperty<MatrixWorkspace>("SecondTransmissionRun", "", Direction::Input,
PropertyMode::Optional, inputValidator->clone()),
"Second, high wavelength transmission run. Optional. Causes the FirstTransmissionRun to be treated as the low wavelength transmission run.");

this->initStitchingInputs();

declareProperty(
new PropertyWithValue<bool>("StrictSpectrumChecking", true, Direction::Input),
"Enforces spectrum number checking prior to normalisation");
declareProperty(new PropertyWithValue<bool>("StrictSpectrumChecking", true, Direction::Input),
"Enforces spectrum number checking prior to normalisation");

setPropertyGroup("FirstTransmissionRun", "Transmission");
setPropertyGroup("SecondTransmissionRun", "Transmission");
Expand All @@ -212,47 +230,83 @@ namespace Mantid

/**
* Correct the position of the detectors based on the input theta value.
* @param toCorrect : Workspace to correct detector posisitions on.
* @param toCorrect : Workspace to correct detector positions on.
* @param thetaInDeg : Theta in degrees to use in correction calculations.
* @param sample : Pointer to the sample
* @param detector : Pointer to a given detector
* @param isPointDetector : True if using point detector analysis
*/
void ReflectometryReductionOne::correctPosition(API::MatrixWorkspace_sptr toCorrect,
const double& thetaInDeg, IComponent_const_sptr sample, IComponent_const_sptr detector)
const double& thetaInDeg, IComponent_const_sptr sample, IComponent_const_sptr detector,
const bool isPointDetector)
{
/*
auto instrument = toCorrect->getInstrument();
const V3D detectorPosition = detector->getPos();
const V3D samplePosition = sample->getPos();
const V3D sampleToDetector = detectorPosition - samplePosition;
auto instrument = toCorrect->getInstrument();
auto referenceFrame = instrument->getReferenceFrame();
const V3D detectorPosition = detector->getPos();
const double sampleToDetectorAlongBeam = sampleToDetector.scalar_prod(
referenceFrame->vecPointingAlongBeam());
const V3D samplePosition = sample->getPos();
const double thetaInRad = thetaInDeg * (M_PI / 180.0);
const V3D sampleToDetector = detectorPosition - samplePosition;
double acrossOffset = 0;
auto referenceFrame = instrument->getReferenceFrame();
double beamOffset = detectorPosition.scalar_prod(referenceFrame->vecPointingAlongBeam());
const double sampleToDetectorAlongBeam = sampleToDetector.scalar_prod(
referenceFrame->vecPointingAlongBeam());
double upOffset = sampleToDetectorAlongBeam * std::sin(2.0 * thetaInRad);
const double thetaInRad = thetaInDeg * (M_PI / 180.0);
auto moveComponentAlg = this->createChildAlgorithm("MoveInstrumentComponent");
moveComponentAlg->initialize();
moveComponentAlg->setProperty("Workspace", toCorrect);
moveComponentAlg->setProperty("ComponentName", detector->getName());
moveComponentAlg->setProperty("RelativePosition", false);
// Movements
moveComponentAlg->setProperty(referenceFrame->pointingAlongBeamAxis(), beamOffset);
moveComponentAlg->setProperty(referenceFrame->pointingHorizontalAxis(), acrossOffset);
moveComponentAlg->setProperty(referenceFrame->pointingUpAxis(), upOffset);
// Execute the movement.
moveComponentAlg->execute();
*/

g_log.debug("Correcting position using theta.");

double acrossOffset = 0;
auto correctPosAlg = this->createChildAlgorithm("SpecularReflectionPositionCorrect");
correctPosAlg->initialize();
correctPosAlg->setProperty("InputWorkspace", toCorrect);

double beamOffset = detectorPosition.scalar_prod(referenceFrame->vecPointingAlongBeam());
const std::string analysisMode = this->getProperty("AnalysisMode");
correctPosAlg->setProperty("AnalysisMode", analysisMode);
const std::string sampleComponentName = this->getProperty("SampleComponentName");
correctPosAlg->setProperty("SampleComponentName", sample->getName());
correctPosAlg->setProperty("TwoThetaIn", thetaInDeg * 2);

double upOffset = sampleToDetectorAlongBeam * std::sin(2.0 * thetaInRad);
if (isPointDetector)
{
correctPosAlg->setProperty("DetectorComponentName", detector->getName());
}
else
{
auto specNumbers = getSpectrumNumbers(toCorrect);
g_log.notice("In Order?");
std::stringstream buffer;
for(size_t t = 0; t < specNumbers.size();++t)
{
buffer << specNumbers[t] << std::endl;
}
g_log.notice(buffer.str());

auto moveComponentAlg = this->createChildAlgorithm("MoveInstrumentComponent");
moveComponentAlg->initialize();
moveComponentAlg->setProperty("Workspace", toCorrect);
moveComponentAlg->setProperty("ComponentName", detector->getName());
moveComponentAlg->setProperty("RelativePosition", false);
// Movements
moveComponentAlg->setProperty(referenceFrame->pointingAlongBeamAxis(), beamOffset);
moveComponentAlg->setProperty(referenceFrame->pointingHorizontalAxis(), acrossOffset);
moveComponentAlg->setProperty(referenceFrame->pointingUpAxis(), upOffset);
// Execute the movement.
moveComponentAlg->execute();
correctPosAlg->setProperty("SpectrumNumbersOfDetectors", specNumbers);
}
correctPosAlg->execute();
MatrixWorkspace_sptr outWS = correctPosAlg->getProperty("OutputWorkspace");
toCorrect = outWS;

}

Expand All @@ -266,7 +320,8 @@ namespace Mantid
* @return
*/
Mantid::API::MatrixWorkspace_sptr ReflectometryReductionOne::toIvsQ(
API::MatrixWorkspace_sptr toConvert, const bool bCorrectPosition, OptionalDouble& thetaInDeg, const bool isPointDetector)
API::MatrixWorkspace_sptr toConvert, const bool bCorrectPosition, OptionalDouble& thetaInDeg,
const bool isPointDetector)
{
/*
* Can either calculate a missing theta value for the purposes of reporting, or correct positions based on a theta value,
Expand All @@ -279,16 +334,24 @@ namespace Mantid
auto correctThetaAlg = this->createChildAlgorithm("SpecularReflectionCalculateTheta");
correctThetaAlg->initialize();
correctThetaAlg->setProperty("InputWorkspace", toConvert);
const std::string detectorComponentName = this->getPropertyValue("DetectorComponentName");
const std::string sampleComponentName = this->getProperty("SampleComponentName");
const std::string analysisMode = this->getProperty("AnalysisMode");
correctThetaAlg->setProperty("DetectorComponentName", detectorComponentName);
correctThetaAlg->setProperty("SampleComponentName", sampleComponentName);
correctThetaAlg->setProperty("AnalysisMode", analysisMode);
const std::string sampleComponentName = this->getProperty("SampleComponentName");
correctThetaAlg->setProperty("SampleComponentName", sampleComponentName);
if (isPointDetector)
{
const std::string detectorComponentName = this->getPropertyValue("DetectorComponentName");
correctThetaAlg->setProperty("DetectorComponentName", detectorComponentName);
}
else
{
std::vector<int> spectrumNumbers = getSpectrumNumbers(toConvert);
correctThetaAlg->setProperty("SpectrumNumbersOfDetectors", spectrumNumbers);
}
correctThetaAlg->execute();
const double twoTheta = correctThetaAlg->getProperty("TwoTheta");

thetaInDeg = twoTheta/2;
thetaInDeg = twoTheta / 2;

}
else if (bCorrectPosition) // This probably ought to be an automatic decision. How about making a guess about sample position holder and detector names. But also allowing the two component names (sample and detector) to be passed in.
Expand All @@ -299,7 +362,7 @@ namespace Mantid
IComponent_const_sptr detector = this->getDetectorComponent(instrument, isPointDetector);
IComponent_const_sptr sample = this->getSurfaceSampleComponent(instrument);

correctPosition(toConvert, thetaInDeg.get(), sample, detector);
correctPosition(toConvert, thetaInDeg.get(), sample, detector, isPointDetector);
}

// Always convert units.
Expand Down Expand Up @@ -471,7 +534,8 @@ namespace Mantid
IvsLam = this->transmissonCorrection(IvsLam, wavelengthInterval,
monitorBackgroundWavelengthInterval, monitorIntegrationWavelengthInterval, i0MonitorIndex,
firstTransmissionRun.get(), secondTransmissionRun, stitchingStart, stitchingDelta,
stitchingEnd, stitchingStartOverlap, stitchingEndOverlap, wavelengthStep, processingCommands);
stitchingEnd, stitchingStartOverlap, stitchingEndOverlap, wavelengthStep,
processingCommands);
}
else
{
Expand Down Expand Up @@ -510,25 +574,26 @@ namespace Mantid
MatrixWorkspace_sptr firstTransmissionRun, OptionalMatrixWorkspace_sptr secondTransmissionRun,
const OptionalDouble& stitchingStart, const OptionalDouble& stitchingDelta,
const OptionalDouble& stitchingEnd, const OptionalDouble& stitchingStartOverlap,
const OptionalDouble& stitchingEndOverlap, const double& wavelengthStep, const std::string& numeratorProcessingCommands)
const OptionalDouble& stitchingEndOverlap, const double& wavelengthStep,
const std::string& numeratorProcessingCommands)
{
g_log.debug("Extracting first transmission run workspace indexes from spectra");

const bool strictSpectrumChecking = getProperty("StrictSpectrumChecking");

MatrixWorkspace_sptr denominator = firstTransmissionRun;
Unit_const_sptr xUnit = firstTransmissionRun->getAxis(0)->unit();
if (xUnit->unitID() == tofUnitId)
{
std::string spectrumProcessingCommands = numeratorProcessingCommands;
/*
If we have strict spectrum checking, the processing commands need to be made from the
numerator workspace AND the transmission workspace based on matching spectrum numbers.
*/
if(strictSpectrumChecking)
If we have strict spectrum checking, the processing commands need to be made from the
numerator workspace AND the transmission workspace based on matching spectrum numbers.
*/
if (strictSpectrumChecking)
{
spectrumProcessingCommands = createWorkspaceIndexListFromDetectorWorkspace(IvsLam,
firstTransmissionRun);
firstTransmissionRun);
}

// Make the transmission run.
Expand Down Expand Up @@ -575,11 +640,12 @@ namespace Mantid
}

/**
@param ws1 : First workspace to compare
@param ws2 : Second workspace to compare against
@param severe: True to indicate that failure to verify should result in an exception. Otherwise a warning is generated.
*/
void ReflectometryReductionOne::verifySpectrumMaps(MatrixWorkspace_const_sptr ws1, MatrixWorkspace_const_sptr ws2, const bool severe)
@param ws1 : First workspace to compare
@param ws2 : Second workspace to compare against
@param severe: True to indicate that failure to verify should result in an exception. Otherwise a warning is generated.
*/
void ReflectometryReductionOne::verifySpectrumMaps(MatrixWorkspace_const_sptr ws1,
MatrixWorkspace_const_sptr ws2, const bool severe)
{
auto map1 = ws1->getSpectrumToWorkspaceIndexMap();
auto map2 = ws2->getSpectrumToWorkspaceIndexMap();
Expand Down
24 changes: 22 additions & 2 deletions ...ework/PythonInterface/test/python/plugins/algorithms/ReflectometryReductionOneBaseTest.py
View file
Expand Up @@ -36,8 +36,7 @@ def setUp(self):
def tearDown(self):
DeleteWorkspace(self.__tof)
DeleteWorkspace(self.__not_tof)



def test_check_input_workpace_not_tof_throws(self):
alg = self.construct_standard_algorithm()
alg.set_InputWorkspace(self.__not_tof)
Expand Down Expand Up @@ -262,12 +261,33 @@ def test_calculate_theta(self):
self.assertAlmostEqual(0.70969419, theta, 4)

DeleteWorkspace(real_run)

def test_correct_positions_point_detector(self):
alg = self.construct_standard_algorithm()
real_run = Load('INTER00013460.nxs')
alg.set_InputWorkspace(real_run)
alg.set_ProcessingInstructions("3,4")
alg.set_FirstTransmissionRun(real_run) # Currently a requirement that one transmisson correction is provided.
alg.set_ThetaIn(0.4) # Low angle
alg.set_CorrectDetectorPositions(True)
out_ws_q1, out_ws_lam1, theta1 = alg.execute()
pos1 = out_ws_lam1.getInstrument().getComponentByName('point-detector').getPos()

alg.set_ThetaIn(0.8) # Repeat with greater incident angle
out_ws_q2, out_ws_lam2, theta2 = alg.execute()
pos2 = out_ws_lam2.getInstrument().getComponentByName('point-detector').getPos()

self.assertTrue(pos2.Y() > pos1.Y(), "Greater incident angle so greater height.")
self.assertEqual(pos2.X(), pos1.X())
self.assertEqual(pos2.Z(), pos1.Z())
DeleteWorkspace(real_run)

def test_multidetector_run(self):
alg = self.construct_standard_algorithm()
real_run = Load('POLREF00004699.nxs')
alg.set_InputWorkspace(real_run[0])
alg.set_AnalysisMode("MultiDetectorAnalysis")
alg.set_CorrectDetectorPositions(False)
alg.set_DetectorComponentName('lineardetector')
alg.set_ProcessingInstructions("3, 10") # Fictional values
alg.set_RegionOfDirectBeam("20, 30") # Fictional values
Expand Down

0 comments on commit d0d12dd

Please sign in to comment.