From d0d12dd1ffc5e0c102d319b21acfde55f23527b2 Mon Sep 17 00:00:00 2001 From: Owen Arnold Date: Mon, 2 Jun 2014 16:08:45 +0100 Subject: [PATCH] refs #9184. Test point detector position correction. --- .../ReflectometryReductionOne.h | 2 +- .../src/ReflectometryReductionOne.cpp | 194 ++++++++++++------ .../ReflectometryReductionOneBaseTest.py | 24 ++- 3 files changed, 153 insertions(+), 67 deletions(-) diff --git a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ReflectometryReductionOne.h b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ReflectometryReductionOne.h index 3e94ce1ece09..b49e9a979bd7 100644 --- a/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ReflectometryReductionOne.h +++ b/Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/ReflectometryReductionOne.h @@ -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); diff --git a/Code/Mantid/Framework/Algorithms/src/ReflectometryReductionOne.cpp b/Code/Mantid/Framework/Algorithms/src/ReflectometryReductionOne.cpp index 869fb37b30be..9fc1da1a60cc 100644 --- a/Code/Mantid/Framework/Algorithms/src/ReflectometryReductionOne.cpp +++ b/Code/Mantid/Framework/Algorithms/src/ReflectometryReductionOne.cpp @@ -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*/ @@ -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 getSpectrumNumbers(MatrixWorkspace_sptr& ws) + { + auto specToWSIndexMap = ws->getSpectrumToWorkspaceIndexMap(); + std::vector keys(specToWSIndexMap.size()); + size_t i = 0; + for (auto it = specToWSIndexMap.begin(); it != specToWSIndexMap.end(); ++it, ++i) + { + keys[i] = static_cast(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 */ @@ -175,7 +194,7 @@ namespace Mantid declareProperty( new WorkspaceProperty("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("SecondTransmissionRun", "", Direction::Input, PropertyMode::Optional, inputValidator->clone()), @@ -183,9 +202,8 @@ namespace Mantid this->initStitchingInputs(); - declareProperty( - new PropertyWithValue("StrictSpectrumChecking", true, Direction::Input), - "Enforces spectrum number checking prior to normalisation"); + declareProperty(new PropertyWithValue("StrictSpectrumChecking", true, Direction::Input), + "Enforces spectrum number checking prior to normalisation"); setPropertyGroup("FirstTransmissionRun", "Transmission"); setPropertyGroup("SecondTransmissionRun", "Transmission"); @@ -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; } @@ -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, @@ -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 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. @@ -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. @@ -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 { @@ -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. @@ -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(); diff --git a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/ReflectometryReductionOneBaseTest.py b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/ReflectometryReductionOneBaseTest.py index 9b7c6e6135f8..3330703297bc 100644 --- a/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/ReflectometryReductionOneBaseTest.py +++ b/Code/Mantid/Framework/PythonInterface/test/python/plugins/algorithms/ReflectometryReductionOneBaseTest.py @@ -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) @@ -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