diff --git a/Framework/API/src/WorkspaceFactory.cpp b/Framework/API/src/WorkspaceFactory.cpp
index e2f890671126a479d732b36c980a5f9a530496d7..2eaca592865c4e520599f043bdb02b4c632c5acd 100644
--- a/Framework/API/src/WorkspaceFactory.cpp
+++ b/Framework/API/src/WorkspaceFactory.cpp
@@ -86,7 +86,7 @@ WorkspaceFactoryImpl::create(const MatrixWorkspace_const_sptr &parent,
return ws;
}
-/** Initialise a workspace from its parent
+/** Initialize a workspace from its parent
* This sets values such as title, instrument, units, sample, spectramap.
* This does NOT copy any data.
*
diff --git a/Framework/Algorithms/CMakeLists.txt b/Framework/Algorithms/CMakeLists.txt
index 445f1c602b5a55ebdbf14fe468ff005ace9d28bc..9231a99e1f936d0aeae54af68af3dcbb6cf7b0f6 100644
--- a/Framework/Algorithms/CMakeLists.txt
+++ b/Framework/Algorithms/CMakeLists.txt
@@ -128,8 +128,10 @@ set ( SRC_FILES
src/GenerateEventsFilter.cpp
src/GeneratePeaks.cpp
src/GeneratePythonScript.cpp
+ src/GetAllEi.cpp
src/GetDetOffsetsMultiPeaks.cpp
src/GetDetectorOffsets.cpp
+ src/GetAllEi.cpp
src/GetEi.cpp
src/GetEi2.cpp
src/GetTimeSeriesLogInformation.cpp
@@ -399,6 +401,7 @@ set ( INC_FILES
inc/MantidAlgorithms/FixGSASInstrumentFile.h
inc/MantidAlgorithms/FlatPlateAbsorption.h
inc/MantidAlgorithms/GSLFunctions.h
+ inc/MantidAlgorithms/GetAllEi.h
inc/MantidAlgorithms/GeneralisedSecondDifference.h
inc/MantidAlgorithms/GenerateEventsFilter.h
inc/MantidAlgorithms/GeneratePeaks.h
@@ -684,6 +687,7 @@ set ( TEST_FILES
GenerateIPythonNotebookTest.h
GeneratePeaksTest.h
GeneratePythonScriptTest.h
+ GetAllEiTest.h
GetDetOffsetsMultiPeaksTest.h
GetDetectorOffsetsTest.h
GetEiTest.h
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h b/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h
index 4a74b3b6a95bf4bd3076f123cf20baf15b8ff2ee..cef126ba1c1e1e557da6b0708b2930960bd8d1e4 100644
--- a/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h
+++ b/Framework/Algorithms/inc/MantidAlgorithms/CreateSampleWorkspace.h
@@ -76,6 +76,7 @@ private:
double noiseScale);
void replaceAll(std::string &str, const std::string &from,
const std::string &to);
+ void addChopperParameters(API::MatrixWorkspace_sptr &ws);
/// A pointer to the random number generator
Kernel::PseudoRandomNumberGenerator *m_randGen;
diff --git a/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h b/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h
new file mode 100644
index 0000000000000000000000000000000000000000..155d8fc2d45d2c6bfa0310962dcfe8cef81e6614
--- /dev/null
+++ b/Framework/Algorithms/inc/MantidAlgorithms/GetAllEi.h
@@ -0,0 +1,126 @@
+#ifndef MANTID_ALGORITHMS_GETALLEI_H_
+#define MANTID_ALGORITHMS_GETALLEI_H_
+
+#include "MantidKernel/System.h"
+#include "MantidKernel/cow_ptr.h"
+#include "MantidAPI/Algorithm.h"
+#include "MantidAPI/MatrixWorkspace.h"
+//#include "MantidAPI/IAlgorithm.h"
+
+namespace Mantid {
+
+namespace Algorithms {
+
+/** Estimate all incident energies, used by chopper instrument.
+
+ Copyright © 2008-9 ISIS Rutherford Appleton Laboratory, NScD Oak Ridge
+ National Laboratory & European Spallation Source
+
+ This file is part of Mantid.
+
+ Mantid is free software; you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version.
+
+ Mantid is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+ File change history is stored at: <https://github.com/mantidproject/mantid>.
+ Code Documentation is available at: <http://doxygen.mantidproject.org>
+*/
+class DLLExport GetAllEi : public API::Algorithm {
+public:
+ GetAllEi();
+ virtual ~GetAllEi(){};
+
+ /// Algorithms name for identification. @see Algorithm::name
+ virtual const std::string name() const { return "GetAllEi"; };
+ /// Algorithm's summary for use in the GUI and help. @see Algorithm::summary
+ virtual const std::string summary() const {
+ return "Analyze the chopper logs and the signal registered by the monitors "
+ "to identify energies used as incident energies in an inelastic "
+ "experiment.";
+ }
+ /// Algorithm's version for identification. @see Algorithm::version
+ virtual int version() const { return 1; };
+ /// Algorithm's category for identification. @see Algorithm::category
+ virtual const std::string category() const { return "Inelastic"; };
+ /// Cross-check properties with each other @see IAlgorithm::validateInputs
+ virtual std::map<std::string, std::string> validateInputs();
+
+private:
+ // Implement abstract Algorithm methods
+ void init();
+ void exec();
+ Kernel::Property *getPLogForProperty(const API::MatrixWorkspace_sptr &inputWS,
+ const std::string &name);
+ void setFilterLog(const API::MatrixWorkspace_sptr &inputWS);
+ // former lambda function exposed as not evry compiler support this yet
+ bool peakGuess(const API::MatrixWorkspace_sptr &inputWS, size_t index,
+ double Ei, const std::vector<size_t> &monsRangeMin,
+ const std::vector<size_t> &monsRangeMax, double &peakPos,
+ double &peakHeight, double &peakTwoSigma);
+
+protected: // for testing, private otherwise.
+ // prepare working workspace with appropriate monitor spectra for fitting
+ API::MatrixWorkspace_sptr
+ // prepare matrix workspace to analyze monitor signal
+ buildWorkspaceToFit(const API::MatrixWorkspace_sptr &inputWS,
+ size_t &wsIndex0);
+
+ /**Return average time series log value for the appropriately filtered log*/
+ double getAvrgLogValue(const API::MatrixWorkspace_sptr &inputWS,
+ const std::string &propertyName,
+ std::vector<Kernel::SplittingInterval> &splitter);
+ /**process logs and retrieve chopper speed and chopper delay*/
+ void findChopSpeedAndDelay(const API::MatrixWorkspace_sptr &inputWS,
+ double &chop_speed, double &chop_delay);
+ void findGuessOpeningTimes(const std::pair<double, double> &TOF_range,
+ double ChopDelay, double Period,
+ std::vector<double> &guess_opening_times);
+ /**Get energy of monitor peak if one is present*/
+ bool findMonitorPeak(const API::MatrixWorkspace_sptr &inputWS, double Ei,
+ const std::vector<size_t> &monsRangeMin,
+ const std::vector<size_t> &monsRangeMax, double &energy,
+ double &height, double &width);
+ /**Find indexes of each expected peak intervals */
+ void findBinRanges(const MantidVec &eBins, const MantidVec &signal,
+ const std::vector<double> &guess_energies,
+ double Eresolution, std::vector<size_t> &irangeMin,
+ std::vector<size_t> &irangeMax,
+ std::vector<bool> &guessValid);
+
+ size_t calcDerivativeAndCountZeros(const std::vector<double> &bins,
+ const std::vector<double> &signal,
+ std::vector<double> &deriv,
+ std::vector<double> &zeros);
+
+ /// if true, take derivate of the filter log to identify interval when
+ /// instrument is running.
+ bool m_FilterWithDerivative;
+ /// maximal relative peak width to consider acceptable. Defined by minimal
+ /// instrument resolution
+ /// and does not exceed 0.08
+ double m_min_Eresolution;
+ // set as half max LET resolution at 20mev at 5e-4
+ double m_max_Eresolution;
+ double m_peakEnergyRatio2reject;
+ // the value of constant phase shift on the chopper used to calculate
+ // tof at chopper from recorded delay.
+ double m_phase;
+ // internal pointer to access to chopper
+ boost::shared_ptr<const Geometry::IComponent> m_chopper;
+ // internal pointer to access log, used for filtering
+ Kernel::TimeSeriesProperty<double> *m_pFilterLog;
+};
+
+} // namespace Algorithms
+} // namespace Mantid
+
+#endif /* MANTID_ALGORITHMS_GETALLEI_H_ */
diff --git a/Framework/Algorithms/src/CreateSampleWorkspace.cpp b/Framework/Algorithms/src/CreateSampleWorkspace.cpp
index 8d3585304e8c43cf9d972af12777e7adeb20ae92..4976102fe9062d9729b35fd29b8bb71cf05e25f1 100644
--- a/Framework/Algorithms/src/CreateSampleWorkspace.cpp
+++ b/Framework/Algorithms/src/CreateSampleWorkspace.cpp
@@ -203,7 +203,6 @@ void CreateSampleWorkspace::exec() {
}
m_randGen = new Kernel::MersenneTwister(seedValue);
}
-
// Create an instrument with one or more rectangular banks.
Instrument_sptr inst = createTestInstrumentRectangular(
numBanks, bankPixelWidth, pixelSpacing, bankDistanceFromSample,
@@ -221,6 +220,8 @@ void CreateSampleWorkspace::exec() {
numBanks * bankPixelWidth * bankPixelWidth, num_bins, xMin, binWidth,
bankPixelWidth * bankPixelWidth, inst, functionString, isRandom);
}
+ // add chopper
+ this->addChopperParameters(ws);
// Set the Unit of the X Axis
try {
@@ -249,6 +250,34 @@ void CreateSampleWorkspace::exec() {
setProperty("OutputWorkspace", ws);
;
}
+/** Add chopper to the existing matrix workspace
+@param ws -- shared pointer to existing matrix workspace which has instrument
+and chopper
+
+@returns workspace modified to have Fermi chopper added to it.
+*/
+void CreateSampleWorkspace::addChopperParameters(
+ API::MatrixWorkspace_sptr &ws) {
+
+ auto testInst = ws->getInstrument();
+ auto chopper = testInst->getComponentByName("chopper-position");
+
+ // add chopper parameters
+ auto ¶mMap = ws->instrumentParameters();
+ const std::string description(
+ "The initial rotation phase of the disk used to calculate the time"
+ " for neutrons arriving at the chopper according to the formula time = "
+ "delay + initial_phase/Speed");
+ paramMap.add<double>("double", chopper.get(), "initial_phase", -3000.,
+ &description);
+ paramMap.add<std::string>("string", chopper.get(), "ChopperDelayLog",
+ "fermi_delay");
+ paramMap.add<std::string>("string", chopper.get(), "ChopperSpeedLog",
+ "fermi_speed");
+ paramMap.add<std::string>("string", chopper.get(), "FilterBaseLog",
+ "is_running");
+ paramMap.add<bool>("bool", chopper.get(), "filter_with_derivative", false);
+}
/** Create histogram workspace
*/
@@ -479,6 +508,12 @@ Instrument_sptr CreateSampleWorkspace::createTestInstrumentRectangular(
testInst->add(source);
testInst->markAsSource(source);
+ // Add chopper
+ ObjComponent *chopper = new ObjComponent(
+ "chopper-position", Object_sptr(new Object), testInst.get());
+ chopper->setPos(V3D(0.0, 0.0, -0.25 * sourceSampleDistance));
+ testInst->add(chopper);
+
// Define a sample as a simple sphere
Object_sptr sampleSphere =
createSphere(0.001, V3D(0.0, 0.0, 0.0), "sample-shape");
diff --git a/Framework/Algorithms/src/FitPeak.cpp b/Framework/Algorithms/src/FitPeak.cpp
index 5a8f34e8eb8dc6a096bccab2c3789419c5e8bac7..831efbb6eeee6d3c66d2fc978891c605bb54ae40 100644
--- a/Framework/Algorithms/src/FitPeak.cpp
+++ b/Framework/Algorithms/src/FitPeak.cpp
@@ -502,10 +502,19 @@ void FitOneSinglePeak::highBkgdFit() {
g_log.warning(outss.str());
size_t numpts = i_maxFitX - i_minFitX;
- i_minPeakX += static_cast<size_t>(static_cast<double>(numpts) / 6.);
- m_minPeakX = m_dataWS->readX(m_wsIndex)[i_minPeakX];
- i_maxPeakX -= static_cast<size_t>(static_cast<double>(numpts) / 6.);
- m_maxPeakX = m_dataWS->readX(m_wsIndex)[i_maxPeakX];
+ size_t shift = static_cast<size_t>(static_cast<double>(numpts) / 6.);
+ i_minPeakX += shift;
+ auto Xdata = m_dataWS->readX(m_wsIndex);
+ if (i_minPeakX >= Xdata.size())
+ i_minPeakX = Xdata.size() - 1;
+ m_minPeakX = Xdata[i_minPeakX];
+
+ if (i_maxPeakX < shift) {
+ i_maxPeakX = 0;
+ } else {
+ i_maxPeakX -= shift;
+ }
+ m_maxPeakX = Xdata[i_maxPeakX];
}
m_bkgdFunc = fitBackground(m_bkgdFunc);
@@ -1288,7 +1297,8 @@ void FitPeak::processProperties() {
// Peak range
vector<double> peakrange = getProperty("PeakRange");
if (peakrange.size() != 2) {
- throw runtime_error("Must enter 2 and only 2 items in fit window. ");
+ throw runtime_error(
+ "Must enter 2 and only 2 items for PeakRange in fit window. ");
}
m_minPeakX = peakrange[0];
m_maxPeakX = peakrange[1];
diff --git a/Framework/Algorithms/src/GetAllEi.cpp b/Framework/Algorithms/src/GetAllEi.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..60418d5071ea437e50237e3f826e2d099efadc9e
--- /dev/null
+++ b/Framework/Algorithms/src/GetAllEi.cpp
@@ -0,0 +1,1255 @@
+//----------------------------------------------------------------------
+// Includes
+//----------------------------------------------------------------------
+#include <boost/format.hpp>
+#include <boost/algorithm/string.hpp>
+#include <string>
+
+#include "MantidAlgorithms/GetAllEi.h"
+#include "MantidKernel/BoundedValidator.h"
+#include "MantidKernel/FilteredTimeSeriesProperty.h"
+#include "MantidKernel/EnabledWhenProperty.h"
+#include "MantidKernel/UnitFactory.h"
+#include "MantidKernel/Unit.h"
+#include "MantidKernel/VectorHelper.h"
+#include "MantidDataObjects/TableWorkspace.h"
+
+namespace Mantid {
+namespace Algorithms {
+
+DECLARE_ALGORITHM(GetAllEi)
+
+/// Empty default constructor
+GetAllEi::GetAllEi()
+ : Algorithm(), m_FilterWithDerivative(true),
+ // minimal resolution for all instruments
+ m_min_Eresolution(0.08),
+ // half maximal resolution for LET
+ m_max_Eresolution(0.5e-3), m_peakEnergyRatio2reject(0.1), m_phase(0),
+ m_chopper(), m_pFilterLog(NULL) {}
+
+/// Initialization method.
+void GetAllEi::init() {
+
+ declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>(
+ "Workspace", "", Kernel::Direction::Input),
+ "The input workspace containing the monitor's spectra "
+ "measured after the last chopper");
+ auto nonNegative = boost::make_shared<Kernel::BoundedValidator<int>>();
+ nonNegative->setLower(0);
+
+ declareProperty(
+ "Monitor1SpecID", EMPTY_INT(), nonNegative,
+ "The workspace index (ID) of the spectra, containing first monitor's"
+ " signal to analyze.");
+ declareProperty(
+ "Monitor2SpecID", EMPTY_INT(), nonNegative,
+ "The workspace index (ID) of the spectra, containing second monitor's"
+ " signal to analyze.");
+
+ declareProperty("ChopperSpeedLog", "Defined in IDF",
+ "Name of the instrument log, "
+ "containing chopper angular velocity. If 'Defined in IDF' "
+ "option is specified, "
+ "the log name is obtained from the IDF");
+ declareProperty(
+ "ChopperDelayLog", "Defined in IDF",
+ "Name of the instrument log, "
+ "containing chopper delay time or chopper phase v.r.t. the pulse time. "
+ "If 'Defined in IDF' option is specified, "
+ "the log name is obtained from IDF");
+ declareProperty(
+ "FilterBaseLog", "Defined in IDF",
+ "Name of the instrument log, "
+ "with positive values indicating that instrument is running\n "
+ "and 0 or negative that it is not.\n"
+ "The log is used to identify time interval to evaluate"
+ " chopper speed and chopper delay which matter.\n"
+ "If such log is not present, average log values are calculated "
+ "within experiment start&end time range.");
+ declareProperty(
+ "FilterWithDerivative", true,
+ "Use derivative of 'FilterBaseLog' "
+ "rather then log values itself to filter invalid time intervals.\n"
+ "Invalid values are then the "
+ "values where the derivative of the log turns zero.\n"
+ "E.g. the 'proton_chage' log grows for each frame "
+ "when instrument is counting and is constant otherwise.");
+ setPropertySettings(
+ "FilterWithDerivative",
+ new Kernel::EnabledWhenProperty("FilterBaseLog",
+ Kernel::ePropertyCriterion::IS_EQUAL_TO,
+ "Defined in IDF"));
+
+ auto maxInRange = boost::make_shared<Kernel::BoundedValidator<double>>();
+ maxInRange->setLower(1.e-6);
+ maxInRange->setUpper(0.1);
+
+ declareProperty("MaxInstrResolution", 0.0005, maxInRange,
+ "The maximal energy resolution possible for an "
+ "instrument at working energies (full width at half "
+ "maximum). \nPeaks, sharper then "
+ "this width are rejected. Accepted limits are: 1e^(-6)-0.1");
+
+ auto minInRange = boost::make_shared<Kernel::BoundedValidator<double>>();
+ minInRange->setLower(0.001);
+ minInRange->setUpper(0.5);
+ declareProperty(
+ "MinInstrResolution", 0.08, minInRange,
+ "The minimal energy resolution possible for an "
+ "instrument at working energies (full width at half maximum).\n"
+ "Peaks broader then this width are rejected. Accepted limits are: "
+ "0.001-0.5");
+
+ auto peakInRange = boost::make_shared<Kernel::BoundedValidator<double>>();
+ peakInRange->setLower(0.0);
+ minInRange->setUpper(1.);
+ declareProperty(
+ "PeaksRatioToReject", 0.1, peakInRange,
+ "Ratio of a peak energy to the maximal energy among all peaks. "
+ "If the ratio is lower then the value specified here, "
+ "peak is treated as insignificant and rejected.\n"
+ "Accepted limits are:0.0 (All accepted) to 1 -- only one peak \n"
+ "(or peaks with max and equal intensity) are accepted.");
+ declareProperty(
+ "IgnoreSecondMonitor", false,
+ "Usually peaks are analyzed and accepted "
+ "only if identified on both monitors. If this property is set to true, "
+ "only first monitor peaks are analyzed.\n"
+ "This is debugging option as getEi has to use both monitors.");
+
+ declareProperty(
+ new API::WorkspaceProperty<API::Workspace>("OutputWorkspace", "",
+ Kernel::Direction::Output),
+ "Name of the output matrix workspace, containing single spectra with"
+ " monitor peaks energies\n"
+ "together with total intensity within each peak.");
+}
+
+// unnamed namespace for auxiliary file-based compilation units
+namespace {
+
+/**Simple template function to remove invalid data from vector
+*@param guessValid -- boolean vector of indicating if each particular guess is
+*valid
+*@param guess -- vector guess values at input and values with removing
+* invalid parameters at output
+*/
+template <class T>
+void removeInvalidValues(const std::vector<bool> &guessValid,
+ std::vector<T> &guess) {
+ std::vector<T> new_guess;
+ new_guess.reserve(guess.size());
+
+ for (size_t i = 0; i < guessValid.size(); i++) {
+ if (guessValid[i]) {
+ new_guess.push_back(guess[i]);
+ }
+ }
+ new_guess.swap(guess);
+}
+/**Internal class to contain peak information */
+struct peakKeeper {
+ double position;
+ double height;
+ double sigma;
+ double energy;
+
+ peakKeeper(double pos, double heigh, double sig)
+ : position(pos), height(heigh), sigma(sig) {
+ this->energy = std::sqrt(2 * M_PI) * height * sigma;
+ }
+ // to sort peaks
+ bool operator<(const peakKeeper &str) const { return (energy > str.energy); }
+};
+
+} // END unnamed namespace for auxiliary file-based compilation units
+
+/** Executes the algorithm -- found all existing monitor peaks. */
+void GetAllEi::exec() {
+ // Get pointers to the workspace, parameter map and table
+ API::MatrixWorkspace_sptr inputWS = getProperty("Workspace");
+ m_min_Eresolution = getProperty("MinInstrResolution");
+ m_max_Eresolution = getProperty("MaxInstrResolution");
+ m_peakEnergyRatio2reject = getProperty("PeaksRatioToReject");
+
+ ////---> recalculate chopper delay to monitor position:
+ auto pInstrument = inputWS->getInstrument();
+ // auto lastChopPositionComponent =
+ // pInstrument->getComponentByName("chopper-position");
+ // auto chopPoint1 = pInstrument->getChopperPoint(0); ->TODO: BUG! this
+ // operation loses parameters map.
+ m_chopper = pInstrument->getComponentByName("chopper-position");
+ if (!m_chopper)
+ throw std::runtime_error("Instrument " + pInstrument->getName() +
+ " does not have 'chopper-position' component");
+
+ auto phase = m_chopper->getNumberParameter("initial_phase");
+
+ if (phase.size() == 0) {
+ throw std::runtime_error("Can not find initial_phase parameter"
+ " attached to the chopper-position component");
+ }
+ if (phase.size() > 1) {
+ throw std::runtime_error(
+ "Can not deal with multiple phases for initial_phase"
+ " parameter attached to the chopper-position component");
+ }
+ m_phase = phase[0];
+
+ this->setFilterLog(inputWS);
+
+ // auto chopPoint1 = pInstrument->getComponentByName("fermi-chopper");
+ // auto par = chopPoint1->getDoubleParameter("Delay (us)");
+ double chopSpeed, chopDelay;
+ findChopSpeedAndDelay(inputWS, chopSpeed, chopDelay);
+ g_log.debug() << boost::str(
+ boost::format("*Identified avrg ChopSpeed: %8.2f and Delay: %8.2f\n") %
+ chopSpeed % chopDelay);
+
+ auto moderator = pInstrument->getSource();
+ double chopDistance = m_chopper->getDistance(
+ *moderator); // location[0].distance(moderator->getPos());
+ double velocity = chopDistance / chopDelay;
+
+ // build workspace to find monitor's peaks
+ size_t det1WSIndex;
+ auto monitorWS = buildWorkspaceToFit(inputWS, det1WSIndex);
+
+ // recalculate delay time from chopper position to monitor position
+ auto detector1 = inputWS->getDetector(det1WSIndex);
+ double mon1Distance = detector1->getDistance(*moderator);
+ double TOF0 = mon1Distance / velocity;
+
+ //--->> below is reserved until full chopper's implementation is available;
+ // auto nChoppers = pInstrument->getNumberOfChopperPoints();
+ // get last chopper.
+ /*
+ if( nChoppers==0)throw std::runtime_error("Instrument does not have any
+ choppers defined");
+
+ auto lastChopper = pInstrument->getChopperPoint(nChoppers-1);
+ ///<---------------------------------------------------
+ */
+ auto baseSpectrum = inputWS->getSpectrum(det1WSIndex);
+ std::pair<double, double> TOF_range = baseSpectrum->getXDataRange();
+
+ double Period =
+ (0.5 * 1.e+6) / chopSpeed; // 0.5 because some choppers open twice.
+ // Would be nice to have it 1 or 0.5 depending on chopper type, but
+ // it looks like not enough information on what chopper is available on ws;
+ std::vector<double> guess_opening;
+ this->findGuessOpeningTimes(TOF_range, TOF0, Period, guess_opening);
+ if (guess_opening.size() == 0) {
+ throw std::runtime_error(
+ "Can not find any chopper opening time within TOF range: " +
+ boost::lexical_cast<std::string>(TOF_range.first) + ':' +
+ boost::lexical_cast<std::string>(TOF_range.second));
+ } else {
+ g_log.debug() << "*Found : " << guess_opening.size()
+ << " chopper prospective opening within time frame: "
+ << TOF_range.first << " to: " << TOF_range.second
+ << std::endl;
+ g_log.debug() << " Timings are:\n";
+ for (size_t i = 0; i < guess_opening.size(); i++) {
+ g_log.debug() << boost::str(boost::format(" %8.2f; ") % guess_opening[i]);
+ }
+ g_log.debug() << std::endl;
+ }
+ std::pair<double, double> Mon1_Erange =
+ monitorWS->getSpectrum(0)->getXDataRange();
+ std::pair<double, double> Mon2_Erange =
+ monitorWS->getSpectrum(1)->getXDataRange();
+ double eMin = std::max(Mon1_Erange.first, Mon2_Erange.first);
+ double eMax = std::min(Mon1_Erange.second, Mon2_Erange.second);
+ g_log.debug() << boost::str(
+ boost::format(
+ "Monitors record data in energy range Emin=%8.2f; Emax=%8.2f\n") %
+ eMin % eMax);
+
+ // convert to energy
+ std::vector<double> guess_ei;
+ guess_ei.reserve(guess_opening.size());
+ double unused(0.0);
+ auto destUnit = Kernel::UnitFactory::Instance().create("Energy");
+ destUnit->initialize(mon1Distance, 0., 0.,
+ static_cast<int>(Kernel::DeltaEMode::Elastic), 0.,
+ unused);
+ for (size_t i = 0; i < guess_opening.size(); i++) {
+ double eGuess = destUnit->singleFromTOF(guess_opening[i]);
+ if (eGuess > eMin && eGuess < eMax) {
+ guess_ei.push_back(eGuess);
+ }
+ }
+ g_log.debug() << "*From all chopper opening only: " +
+ boost::lexical_cast<std::string>(guess_ei.size()) +
+ " fell within both monitor's recording energy range\n";
+ g_log.debug() << " Guess Energies are:\n";
+ for (size_t i = 0; i < guess_ei.size(); i++) {
+ g_log.debug() << boost::str(boost::format(" %8.2f; ") % guess_ei[i]);
+ }
+ g_log.debug() << std::endl;
+
+ std::sort(guess_ei.begin(), guess_ei.end());
+
+ std::vector<size_t> irange_min, irange_max;
+ std::vector<bool> guessValid;
+ // preprocess first monitors peaks;
+ g_log.debug() << "*Looking for real energy peaks on first monitor\n";
+ findBinRanges(monitorWS->readX(0), monitorWS->readY(0), guess_ei,
+ this->m_min_Eresolution / (2 * std::sqrt(2 * std::log(2.))),
+ irange_min, irange_max, guessValid);
+
+ // remove invalid guess values
+ removeInvalidValues<double>(guessValid, guess_ei);
+
+ // preprocess second monitors peaks
+ std::vector<size_t> irange1_min, irange1_max;
+ if (!this->getProperty("IgnoreSecondMonitor")) {
+ g_log.debug() << "*Looking for real energy peaks on second monitor\n";
+ findBinRanges(monitorWS->readX(1), monitorWS->readY(1), guess_ei,
+ this->m_min_Eresolution / (2 * std::sqrt(2 * std::log(2.))),
+ irange1_min, irange1_max, guessValid);
+ removeInvalidValues<double>(guessValid, guess_ei);
+ removeInvalidValues<size_t>(guessValid, irange_min);
+ removeInvalidValues<size_t>(guessValid, irange_max);
+ } else {
+ // this is wrong but will not be used anyway
+ // (except formally looping through vector)
+ irange1_min.assign(irange_min.begin(), irange_min.end());
+ irange1_max.assign(irange_max.begin(), irange_max.end());
+ }
+ g_log.debug()
+ << "*Identified: " + boost::lexical_cast<std::string>(guess_ei.size()) +
+ " peaks with sufficient signal around guess chopper opening\n";
+
+ std::vector<peakKeeper> peaks;
+
+ double maxPeakEnergy(0);
+ std::vector<size_t> monsRangeMin(2), monsRangeMax(2);
+ for (size_t i = 0; i < guess_ei.size(); i++) {
+ monsRangeMin[0] = irange_min[i];
+ monsRangeMax[0] = irange_max[i];
+ monsRangeMin[1] = irange1_min[i];
+ monsRangeMax[1] = irange1_max[i];
+
+ double energy, height, twoSigma;
+ bool found = findMonitorPeak(monitorWS, guess_ei[i], monsRangeMin,
+ monsRangeMax, energy, height, twoSigma);
+ if (found) {
+ peaks.push_back(peakKeeper(energy, height, 0.5 * twoSigma));
+ if (peaks.back().energy > maxPeakEnergy)
+ maxPeakEnergy = peaks.back().energy;
+ }
+ }
+ monitorWS.reset();
+
+ size_t nPeaks = peaks.size();
+ if (nPeaks == 0) {
+ throw std::runtime_error("Can not identify any energy peaks");
+ }
+ // sort peaks and remove invalid one
+ guessValid.resize(nPeaks);
+ bool needsRemoval(false);
+ for (size_t i = 0; i < nPeaks; i++) {
+ peaks[i].energy /= maxPeakEnergy;
+ if (peaks[i].energy < m_peakEnergyRatio2reject) {
+ guessValid[i] = false;
+ g_log.debug() << "*Rejecting peak at Ei=" +
+ boost::lexical_cast<std::string>(peaks[i].position) +
+ " as its total energy lower then the threshold\n";
+ needsRemoval = true;
+ } else {
+ guessValid[i] = true;
+ }
+ }
+ if (needsRemoval)
+ removeInvalidValues<peakKeeper>(guessValid, peaks);
+ nPeaks = peaks.size();
+ // sort by energy decreasing -- see class definition
+ std::sort(peaks.begin(), peaks.end());
+
+ // finalize output
+ auto result_ws = API::WorkspaceFactory::Instance().create("Workspace2D", 1,
+ nPeaks, nPeaks);
+
+ MantidVec peaks_positions;
+ MantidVec &Signal = result_ws->dataY(0);
+ MantidVec &Error = result_ws->dataE(0);
+ for (size_t i = 0; i < nPeaks; i++) {
+ peaks_positions.push_back(peaks[i].position);
+ Signal[i] = peaks[i].height;
+ Error[i] = peaks[i].sigma;
+ }
+ result_ws->setX(0, peaks_positions);
+
+ setProperty("OutputWorkspace", result_ws);
+}
+
+// unnamed namespace for auxiliary file-based functions, converted from lambda
+// as not all Mantid compilers support lambda yet.
+
+/**The internal procedure to set filter log from properties,
+ defining it.
+* @param inputWS -- shared pointer to the input workspace with
+ logs to analyze
+*/
+void GetAllEi::setFilterLog(const API::MatrixWorkspace_sptr &inputWS) {
+
+ std::string filerLogName;
+ std::string filterBase = getProperty("FilterBaseLog");
+ if (boost::iequals(filterBase, "Defined in IDF")) {
+ filerLogName = m_chopper->getStringParameter("FilterBaseLog")[0];
+ m_FilterWithDerivative =
+ m_chopper->getBoolParameter("filter_with_derivative")[0];
+ } else {
+ filerLogName = filterBase;
+ m_FilterWithDerivative = getProperty("FilterWithDerivative");
+ }
+ try {
+ m_pFilterLog = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(
+ inputWS->run().getProperty(filerLogName));
+ } catch (std::runtime_error &) {
+ g_log.warning() << " Can not retrieve (double) filtering log: " +
+ filerLogName +
+ " from current workspace\n"
+ " Using total experiment range to "
+ "find logs averages for chopper parameters\n";
+ m_FilterWithDerivative = false;
+ }
+}
+/**Former lambda to identify guess values for a peak at given index
+* and set up these parameters as input for fitting algorithm
+*
+*@param inputWS -- the workspace to process
+*@param index -- the number of the workspace spectra to process
+*@param Ei -- incident energy
+*@param monsRangeMin -- vector of left boundaries for the subintervals to look
+*for peak
+*@param monsRangeMax -- vector of right boundaries for the subintervals to look
+*for peak
+*
+*@param peakPos -- output energy of the peak
+*@param peakHeight -- output height of the peak assuming Gaussian shape
+*@param peakTwoSigma -- output width of the peak assuming Gaussian shape
+*/
+bool GetAllEi::peakGuess(const API::MatrixWorkspace_sptr &inputWS, size_t index,
+ double Ei, const std::vector<size_t> &monsRangeMin,
+ const std::vector<size_t> &monsRangeMax,
+ double &peakPos, double &peakHeight,
+ double &peakTwoSigma) {
+
+ // calculate sigma from half-width parameters
+ double maxSigma = Ei * m_min_Eresolution / (2 * std::sqrt(2 * std::log(2.)));
+
+ double sMin(std::numeric_limits<double>::max());
+ double sMax(-sMin);
+ double xOfMax(0), dXmax(0);
+ double Intensity(0);
+
+ auto X = inputWS->readX(index);
+ auto S = inputWS->readY(index);
+ size_t ind_min = monsRangeMin[index];
+ size_t ind_max = monsRangeMax[index];
+ // interval too small -- not interested in a peak there
+ if (std::fabs(double(ind_max - ind_min)) < 5)
+ return false;
+
+ // double xMin = X[ind_min];
+ // double xMax = X[ind_max];
+ // size_t ind_Ofmax(ind_min);
+
+ for (size_t i = ind_min; i < ind_max; i++) {
+ double dX = X[i + 1] - X[i];
+ double signal = S[i] / dX;
+ if (signal < sMin)
+ sMin = signal;
+ if (signal > sMax) {
+ sMax = signal;
+ dXmax = dX;
+ xOfMax = X[i];
+ // ind_Ofmax=i;
+ }
+ Intensity += S[i];
+ }
+ // monitor peak should not have just two counts in it.
+ if (sMax * dXmax <= 2)
+ return false;
+ //
+ // size_t SearchAreaSize = ind_max - ind_min;
+
+ double SmoothRange = 2 * maxSigma;
+
+ std::vector<double> SAvrg, binsAvrg;
+ Kernel::VectorHelper::smoothInRange(S, SAvrg, SmoothRange, &X, ind_min,
+ ind_max, &binsAvrg);
+
+ double realPeakPos(xOfMax); // this position is less shifted
+ // due to the skew in averaging formula
+ bool foundRealPeakPos(false);
+ std::vector<double> der1Avrg, der2Avrg, peaks, hillsPos, SAvrg1, binsAvrg1;
+ size_t nPeaks =
+ this->calcDerivativeAndCountZeros(binsAvrg, SAvrg, der1Avrg, peaks);
+ size_t nHills =
+ this->calcDerivativeAndCountZeros(binsAvrg, der1Avrg, der2Avrg, hillsPos);
+ size_t nPrevHills = 2 * nHills;
+ if (nPeaks == 1) {
+ foundRealPeakPos = true;
+ realPeakPos = peaks[0];
+ }
+
+ size_t ic(0), stay_still_count(0);
+ bool iterations_fail(false);
+ while ((nPeaks > 1 || nHills > 2) && (!iterations_fail)) {
+ Kernel::VectorHelper::smoothInRange(SAvrg, SAvrg1, SmoothRange, &binsAvrg,
+ 0, ind_max - ind_min, &binsAvrg1);
+ nPrevHills = nHills;
+
+ nPeaks =
+ this->calcDerivativeAndCountZeros(binsAvrg1, SAvrg1, der1Avrg, peaks);
+ nHills = this->calcDerivativeAndCountZeros(binsAvrg1, der1Avrg, der2Avrg,
+ hillsPos);
+ SAvrg.swap(SAvrg1);
+ binsAvrg.swap(binsAvrg1);
+ if (nPeaks == 1 && !foundRealPeakPos) { // fix first peak position found
+ foundRealPeakPos = true; // as averaging shift peaks on
+ realPeakPos = peaks[0]; // irregular grid.
+ }
+ ic++;
+ if (nPrevHills <= nHills) {
+ stay_still_count++;
+ } else {
+ stay_still_count = 0;
+ }
+ if (ic > 50 || stay_still_count > 3)
+ iterations_fail = true;
+ }
+ if (iterations_fail) {
+ g_log.information() << "*No peak search convergence after " +
+ boost::lexical_cast<std::string>(ic) +
+ " smoothing iterations at still_count: " +
+ boost::lexical_cast<std::string>(
+ stay_still_count) +
+ " Wrong energy or noisy peak at Ei=" +
+ boost::lexical_cast<std::string>(Ei)
+ << std::endl;
+ }
+ g_log.debug() << "*Performed: " + boost::lexical_cast<std::string>(ic) +
+ " averages for spectra " +
+ boost::lexical_cast<std::string>(index) +
+ " at energy: " + boost::lexical_cast<std::string>(Ei) +
+ "\n and found: " +
+ boost::lexical_cast<std::string>(nPeaks) + "peaks and " +
+ boost::lexical_cast<std::string>(nHills) + " hills\n";
+ if (nPeaks != 1) {
+ g_log.debug() << "*Peak rejected as n-peaks !=1 after averaging\n";
+ return false;
+ }
+
+ peakPos = peaks[0];
+ if (nHills > 2) {
+ size_t peakIndex = Kernel::VectorHelper::getBinIndex(hillsPos, peaks[0]);
+ peakTwoSigma = hillsPos[peakIndex + 1] - hillsPos[peakIndex];
+ } else {
+ if (hillsPos.size() == 2) {
+ peakTwoSigma = hillsPos[1] - hillsPos[0];
+ } else {
+ g_log.debug() << "*Peak rejected as averaging gives: " +
+ boost::lexical_cast<std::string>(nPeaks) +
+ " peaks and " +
+ boost::lexical_cast<std::string>(nHills) +
+ " heals\n";
+
+ return false;
+ }
+ }
+ // assuming that averaging conserves intensity and removing linear
+ // background:
+ peakHeight = Intensity / (0.5 * std::sqrt(2. * M_PI) * peakTwoSigma) - sMin;
+ peakPos = realPeakPos;
+
+ return true;
+}
+
+/**Get energy of monitor peak if one is present
+*@param inputWS -- the workspace to process
+*@param Ei -- incident energy
+*@param monsRangeMin -- vector of indexes of left boundaries
+* for the subintervals to look for peak
+*@param monsRangeMax -- vector of indexes of right boundaries
+* for the subintervals to look for peak
+*
+*@param position -- output energy of the peak center.
+*@param height -- output height of the peak assuming Gaussian shape
+*@param twoSigma -- output width of the peak assuming Gaussian shape
+*/
+bool GetAllEi::findMonitorPeak(const API::MatrixWorkspace_sptr &inputWS,
+ double Ei,
+ const std::vector<size_t> &monsRangeMin,
+ const std::vector<size_t> &monsRangeMax,
+ double &position, double &height,
+ double &twoSigma) {
+ // calculate sigma from half-width parameters
+ double maxSigma = Ei * m_min_Eresolution / (2 * std::sqrt(2 * std::log(2.)));
+ double minSigma = Ei * m_max_Eresolution / (2 * std::sqrt(2 * std::log(2.)));
+ //--------------------------------------------------------------------
+ double peak1Pos, peak1TwoSigma, peak1Height;
+ if (!peakGuess(inputWS, 0, Ei, monsRangeMin, monsRangeMax, peak1Pos,
+ peak1Height, peak1TwoSigma))
+ return false;
+ if (0.25 * peak1TwoSigma > maxSigma || peak1TwoSigma < minSigma) {
+ g_log.debug() << "*Rejecting due to width: Peak at mon1 Ei=" +
+ boost::lexical_cast<std::string>(peak1Pos) +
+ " with Height:" +
+ boost::lexical_cast<std::string>(peak1Height) +
+ " and 2*Sigma: " +
+ boost::lexical_cast<std::string>(peak1TwoSigma)
+ << std::endl;
+ return false;
+ }
+
+ if (!this->getProperty("IgnoreSecondMonitor")) {
+ double peak2Pos, peak2TwoSigma, peak2Height;
+ if (!peakGuess(inputWS, 1, Ei, monsRangeMin, monsRangeMax, peak2Pos,
+ peak2Height, peak2TwoSigma))
+ return false;
+ // Let's not check anything except peak position for monitor2, as
+ // its intensity may be very low for some instruments.
+ // if(0.25*peak2TwoSigma>maxSigma||peak2TwoSigma<minSigma)return false;
+
+ // peak in first and second monitors are too far from each other. May be the
+ // instrument
+ // is ill-calibrated but GetEi will probably not find this peak anyway.
+ if (std::fabs(peak1Pos - peak2Pos) >
+ 0.25 * (peak1TwoSigma + peak2TwoSigma)) {
+ g_log.debug()
+ << "*Rejecting due to displacement between Peak at mon1: Ei=" +
+ boost::lexical_cast<std::string>(peak1Pos) + " with Height:" +
+ boost::lexical_cast<std::string>(peak1Height) +
+ " and 2*Sigma: " +
+ boost::lexical_cast<std::string>(peak1TwoSigma) +
+ "\n and Peak at mon2: Ei= " +
+ boost::lexical_cast<std::string>(peak2Pos) + "and height: " +
+ boost::lexical_cast<std::string>(peak1Height) << std::endl;
+
+ return false;
+ }
+ }
+
+ position = peak1Pos;
+ twoSigma = peak1TwoSigma;
+ height = peak1Height;
+
+ return true;
+}
+namespace { // for lambda extracted from calcDerivativeAndCountZeros
+ /**former lambda from calcDerivativeAndCountZeros
+ *estimating if sign have changed from its previous value
+ *@param val -- current function value
+ *@param prevSign -- the sign of the function at previous value
+ */
+bool signChanged(double val, int &prevSign) {
+ int curSign = (val >= 0 ? 1 : -1);
+ bool changed = curSign != prevSign;
+ prevSign = curSign;
+ return changed;
+}
+}
+
+/**Bare-bone function to calculate numerical derivative, and estimate number of
+* zeros this derivative has. The function is assumed to be defined on the
+* the left of a bin range so the derivative is calculated in the same point.
+* No checks are performed for simplicity so data have to be correct
+* form at input.
+*@param bins -- vector of bin boundaries.
+*@param signal -- vector of signal size of bins.size()-1
+*@param deriv -- output vector of numerical derivative
+*@param zeros -- coordinates of found zeros
+
+*@return -- number of zeros, the derivative has in the interval provided.
+*/
+size_t GetAllEi::calcDerivativeAndCountZeros(const std::vector<double> &bins,
+ const std::vector<double> &signal,
+ std::vector<double> &deriv,
+ std::vector<double> &zeros) {
+ size_t nPoints = signal.size();
+ deriv.resize(nPoints);
+ zeros.resize(0);
+
+ std::list<double> funVal;
+ double bin0 = bins[1] - bins[0];
+ double f0 = signal[0] / bin0;
+ double bin1 = bins[2] - bins[1];
+ double f1 = signal[1] / bin1;
+
+ size_t nZeros(0);
+
+ funVal.push_front(f1);
+ deriv[0] = 2 * (f1 - f0) / (bin0 + bin1);
+ int prevSign = (deriv[0] >= 0 ? 1 : -1);
+
+ for (size_t i = 1; i < nPoints - 1; i++) {
+ bin1 = (bins[i + 2] - bins[i + 1]);
+ f1 = signal[i + 1] / bin1;
+ deriv[i] = (f1 - f0) / (bins[i + 2] - bins[i]);
+ f0 = funVal.back();
+ funVal.pop_back();
+ funVal.push_front(f1);
+
+ if (signChanged(deriv[i], prevSign)) {
+ nZeros++;
+ zeros.push_back(0.5 * (bins[i - 1] + bins[i]));
+ }
+ }
+ deriv[nPoints - 1] = 2 * (f1 - f0) / (bin1 + bin0);
+ if (signChanged(deriv[nPoints - 1], prevSign)) {
+ zeros.push_back(bins[nPoints - 1]);
+ nZeros++;
+ }
+
+ return nZeros;
+}
+namespace { // for lambda extracted from findBinRanges
+// get bin range corresponding to the energy range
+void getBinRange(const MantidVec &eBins, double eMin, double eMax,
+ size_t &index_min, size_t &index_max) {
+
+ size_t nBins = eBins.size();
+ if (eMin <= eBins[0]) {
+ index_min = 0;
+ } else {
+ index_min = Kernel::VectorHelper::getBinIndex(eBins, eMin);
+ }
+
+ if (eMax >= eBins[nBins - 1]) {
+ index_max = nBins - 1;
+ } else {
+ index_max = Kernel::VectorHelper::getBinIndex(eBins, eMax) + 1;
+ if (index_max >= nBins)
+ index_max = nBins - 1; // last bin range anyway. Should not happen
+ }
+}
+
+// refine bin range. May need better procedure for this.
+bool refineEGuess(const MantidVec &eBins, const MantidVec &signal,
+ double &eGuess, size_t index_min, size_t index_max) {
+
+ size_t ind_Emax = index_min;
+ double SMax(0);
+ for (size_t i = index_min; i < index_max; i++) {
+ double dX = eBins[i + 1] - eBins[i];
+ double sig = signal[i] / dX;
+ if (sig > SMax) {
+ SMax = sig;
+ ind_Emax = i;
+ }
+ }
+ if (ind_Emax == index_min || ind_Emax == index_max) {
+ return false;
+ }
+ eGuess = 0.5 * (eBins[ind_Emax] + eBins[ind_Emax + 1]);
+ return true;
+}
+
+struct peakKeeper2 {
+ double left_rng;
+ double right_rng;
+ peakKeeper2(){};
+ peakKeeper2(double left, double right) : left_rng(left), right_rng(right) {}
+};
+}
+
+/**Find indexes of each expected peak intervals from monotonous array of ranges.
+*@param eBins -- bin ranges to look through
+*@param signal -- vector of signal in the bins
+*@param guess_energy -- vector of guess energies to look for
+*@param eResolution -- instrument resolution in energy units
+*@param irangeMin -- start indexes of energy intervals in the guess_energies
+* vector.
+*@param irangeMax -- final indexes of energy intervals in the guess_energies
+* vector.
+*@param guessValid -- output boolean vector, which specifies if guess energies
+* in guess_energy vector are valid or not
+*/
+void GetAllEi::findBinRanges(const MantidVec &eBins, const MantidVec &signal,
+ const std::vector<double> &guess_energy,
+ double eResolution, std::vector<size_t> &irangeMin,
+ std::vector<size_t> &irangeMax,
+ std::vector<bool> &guessValid) {
+
+ // size_t nBins = eBins.size();
+ guessValid.resize(guess_energy.size());
+
+ // Do the job
+ size_t ind_min, ind_max;
+ irangeMin.resize(0);
+ irangeMax.resize(0);
+
+ // identify guess bin ranges
+ std::vector<peakKeeper2> guess_peak(guess_energy.size());
+ for (size_t nGuess = 0; nGuess < guess_energy.size(); nGuess++) {
+ double eGuess = guess_energy[nGuess];
+ getBinRange(eBins, eGuess * (1 - 4 * eResolution),
+ eGuess * (1 + 4 * eResolution), ind_min, ind_max);
+ guess_peak[nGuess] = peakKeeper2(eBins[ind_min], eBins[ind_max]);
+ }
+ // verify that the ranges not intercept and refine interceptions
+ for (size_t i = 1; i < guess_energy.size(); i++) {
+ if (guess_peak[i - 1].right_rng > guess_peak[i].left_rng) {
+ double mid_pnt =
+ 0.5 * (guess_peak[i - 1].right_rng + guess_peak[i].left_rng);
+ guess_peak[i - 1].right_rng = mid_pnt;
+ guess_peak[i].left_rng = mid_pnt;
+ }
+ }
+ // identify final bin ranges
+ for (size_t nGuess = 0; nGuess < guess_energy.size(); nGuess++) {
+
+ double eGuess = guess_energy[nGuess];
+ getBinRange(eBins, guess_peak[nGuess].left_rng,
+ guess_peak[nGuess].right_rng, ind_min, ind_max);
+
+ if (refineEGuess(eBins, signal, eGuess, ind_min, ind_max)) {
+ getBinRange(eBins, std::max(guess_peak[nGuess].left_rng,
+ eGuess * (1 - 3 * eResolution)),
+ std::max(guess_peak[nGuess].right_rng,
+ eGuess * (1 + 3 * eResolution)),
+ ind_min, ind_max);
+ irangeMin.push_back(ind_min);
+ irangeMax.push_back(ind_max);
+ guessValid[nGuess] = true;
+ } else {
+ guessValid[nGuess] = false;
+ g_log.debug() << "*Incorrect guess energy: "
+ << boost::lexical_cast<std::string>(eGuess)
+ << " no energy peak found in 4*Sigma range\n";
+ }
+ }
+ // if array decreasing rather then increasing, indexes behave differently.
+ // Will it still work?
+ if (irangeMax.size() > 0) {
+ if (irangeMax[0] < irangeMin[0]) {
+ irangeMax.swap(irangeMin);
+ }
+ }
+}
+
+/**Build 2-spectra workspace in units of energy, used as source
+*to identify actual monitors spectra
+*@param inputWS shared pointer to initial workspace
+*@param wsIndex0 -- returns workspace index for first detector.
+*@return shared pointer to intermediate workspace, in units of energy
+* used to fit monitor's spectra.
+*/
+API::MatrixWorkspace_sptr
+GetAllEi::buildWorkspaceToFit(const API::MatrixWorkspace_sptr &inputWS,
+ size_t &wsIndex0) {
+
+ // at this stage all properties are validated so its safe to access them
+ // without
+ // additional checks.
+ specid_t specNum1 = getProperty("Monitor1SpecID");
+ wsIndex0 = inputWS->getIndexFromSpectrumNumber(specNum1);
+ specid_t specNum2 = getProperty("Monitor2SpecID");
+ size_t wsIndex1 = inputWS->getIndexFromSpectrumNumber(specNum2);
+ auto pSpectr1 = inputWS->getSpectrum(wsIndex0);
+ auto pSpectr2 = inputWS->getSpectrum(wsIndex1);
+ // assuming equally binned ws.
+ // auto bins = inputWS->dataX(wsIndex0);
+ auto bins = pSpectr1->dataX();
+ size_t XLength = bins.size();
+ size_t YLength = inputWS->dataY(wsIndex0).size();
+ auto working_ws =
+ API::WorkspaceFactory::Instance().create(inputWS, 2, XLength, YLength);
+ // copy data --> very bad as implicitly assigns pointer
+ // to bins array and bins array have to exist out of this routine
+ // scope.
+ // This does not matter in this case as below we convert units
+ // which should decouple cow_pointer but very scary operation in
+ // general.
+ working_ws->setX(0, bins);
+ working_ws->setX(1, bins);
+ MantidVec &Signal1 = working_ws->dataY(0);
+ MantidVec &Error1 = working_ws->dataE(0);
+ MantidVec &Signal2 = working_ws->dataY(1);
+ MantidVec &Error2 = working_ws->dataE(1);
+ for (size_t i = 0; i < YLength; i++) {
+ Signal1[i] = inputWS->dataY(wsIndex0)[i];
+ Error1[i] = inputWS->dataE(wsIndex0)[i];
+ Signal2[i] = inputWS->dataY(wsIndex1)[i];
+ Error2[i] = inputWS->dataE(wsIndex1)[i];
+ }
+ // copy detector mapping
+ API::ISpectrum *spectrum = working_ws->getSpectrum(0);
+ spectrum->setSpectrumNo(specNum1);
+ spectrum->clearDetectorIDs();
+ spectrum->addDetectorIDs(pSpectr1->getDetectorIDs());
+ spectrum = working_ws->getSpectrum(1);
+ spectrum->setSpectrumNo(specNum2);
+ spectrum->clearDetectorIDs();
+ spectrum->addDetectorIDs(pSpectr2->getDetectorIDs());
+
+ if (inputWS->getAxis(0)->unit()->caption() != "Energy") {
+ API::IAlgorithm_sptr conv = createChildAlgorithm("ConvertUnits");
+ conv->initialize();
+ conv->setProperty("InputWorkspace", working_ws);
+ conv->setProperty("OutputWorkspace", working_ws);
+ conv->setPropertyValue("Target", "Energy");
+ conv->setPropertyValue("EMode", "Elastic");
+ // conv->setProperty("AlignBins",true); --> throws due to bug in
+ // ConvertUnits
+ conv->execute();
+ }
+
+ return working_ws;
+}
+/**function calculates list of provisional chopper opening times.
+*@param TOF_range -- std::pair containing min and max time, of signal
+* measured on monitors
+*@param ChopDelay -- the time of flight neutrons travel from source
+* to the chopper opening.
+*@param Period -- period of chopper openings
+
+*@param guess_opening_times -- output vector with time values
+* at which neutrons may pass through the chopper.
+*/
+void GetAllEi::findGuessOpeningTimes(const std::pair<double, double> &TOF_range,
+ double ChopDelay, double Period,
+ std::vector<double> &guess_opening_times) {
+
+ if (ChopDelay >= TOF_range.second) {
+ std::string chop = boost::str(boost::format("%.2g") % ChopDelay);
+ std::string t_min = boost::str(boost::format("%.2g") % TOF_range.first);
+ std::string t_max = boost::str(boost::format("%.2g") % TOF_range.second);
+ throw std::runtime_error("Logical error: Chopper opening time: " + chop +
+ " is outside of time interval: " + t_min + ":" +
+ t_max + " of the signal, measured on monitors.");
+ }
+
+ // number of times chopper with specified rotation period opens.
+ size_t n_openings =
+ static_cast<size_t>((TOF_range.second - ChopDelay) / Period) + 1;
+ // number of periods falling outside of the time period, measuring on monitor.
+ size_t n_start(0);
+ if (ChopDelay < TOF_range.first) {
+ n_start = static_cast<size_t>((TOF_range.first - ChopDelay) / Period) + 1;
+ n_openings -= n_start;
+ }
+
+ guess_opening_times.resize(n_openings);
+ for (size_t i = n_start; i < n_openings + n_start; i++) {
+ guess_opening_times[i - n_start] =
+ ChopDelay + static_cast<double>(i) * Period;
+ }
+}
+/**Finds pointer to log value for the property with the name provided
+*
+*@param inputWS -- workspace with logs attached
+*@param propertyName -- name of the property to find log for
+*
+*@return -- pointer to property which contain the log requested or nullptr if
+* no log found or other errors identified. */
+Kernel::Property *
+GetAllEi::getPLogForProperty(const API::MatrixWorkspace_sptr &inputWS,
+ const std::string &propertyName) {
+
+ std::string LogName = this->getProperty(propertyName);
+ if (boost::iequals(LogName, "Defined in IDF")) {
+ auto AllNames = m_chopper->getStringParameter(propertyName);
+ if (AllNames.size() != 1)
+ return NULL;
+ LogName = AllNames[0];
+ }
+ auto pIProperty = (inputWS->run().getProperty(LogName));
+
+ return pIProperty;
+}
+
+/**Return average time series log value for the appropriately filtered log
+* @param inputWS -- shared pointer to the input workspace containing
+* the log to process
+* @param propertyName -- log name
+* @param splitter -- array of Kernel::splittingInterval data, used to
+* filter input events or empty array to use
+* experiment start/end times.
+*/
+double
+GetAllEi::getAvrgLogValue(const API::MatrixWorkspace_sptr &inputWS,
+ const std::string &propertyName,
+ std::vector<Kernel::SplittingInterval> &splitter) {
+
+ auto pIProperty = getPLogForProperty(inputWS, propertyName);
+
+ // this will always provide a defined pointer as this has been verified in
+ // validator.
+ auto pTimeSeries =
+ dynamic_cast<Kernel::TimeSeriesProperty<double> *>(pIProperty);
+
+ if (splitter.size() == 0) {
+ auto TimeStart = inputWS->run().startTime();
+ auto TimeEnd = inputWS->run().endTime();
+ pTimeSeries->filterByTime(TimeStart, TimeEnd);
+ } else {
+ pTimeSeries->filterByTimes(splitter);
+ }
+ if (pTimeSeries->size() == 0) {
+ throw std::runtime_error(
+ "Can not find average value for log defined by property" +
+ propertyName + " As no valid log values are found.");
+ }
+
+ return pTimeSeries->getStatistics().mean;
+}
+namespace { // former lambda function for findChopSpeedAndDelay
+
+/**Select time interval on the basis of previous time interval
+* selection and check if current value gets in the selection
+*
+* @param t_beg -- initial time for current time interval
+* @param t_end -- final time for current time interval
+* @param inSelection -- the boolean indicating if previous interval
+* was selected on input and current selected on
+* output
+* @param startTime -- total selection time start moment
+* @param endTime -- total selection time final moments
+*
+*@return true if selection interval is completed
+* (current interval is not selected) and false otherwise
+*/
+bool SelectInterval(const Kernel::DateAndTime &t_beg,
+ const Kernel::DateAndTime &t_end, double value,
+ bool &inSelection, Kernel::DateAndTime &startTime,
+ Kernel::DateAndTime &endTime) {
+
+ if (value > 0) {
+ if (!inSelection) {
+ startTime = t_beg;
+ }
+ inSelection = true;
+ } else {
+ if (inSelection) {
+ inSelection = false;
+ if (endTime > startTime)
+ return true;
+ }
+ }
+ endTime = t_end;
+ return false;
+}
+}
+/**Analyze chopper logs and identify chopper speed and delay
+@param inputWS -- sp to workspace with attached logs.
+@param chop_speed -- output value for chopper speed in uSec
+@param chop_delay -- output value for chopper delay in uSec
+*/
+void GetAllEi::findChopSpeedAndDelay(const API::MatrixWorkspace_sptr &inputWS,
+ double &chop_speed, double &chop_delay) {
+
+ // TODO: Make it dependent on inputWS time range
+
+ std::vector<Kernel::SplittingInterval> splitter;
+ if (m_pFilterLog) {
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> pDerivative;
+
+ // Define selecting function
+ bool inSelection(false);
+ // time interval to select (start-end)
+ Kernel::DateAndTime startTime, endTime;
+ //
+ // Analyze filtering log
+ auto dateAndTimes = m_pFilterLog->valueAsCorrectMap();
+ auto it = dateAndTimes.begin();
+ auto next = it;
+ next++;
+ std::map<Kernel::DateAndTime, double> derivMap;
+ auto itder = it;
+ if (m_FilterWithDerivative) {
+ pDerivative = m_pFilterLog->getDerivative();
+ derivMap = pDerivative->valueAsCorrectMap();
+ itder = derivMap.begin();
+ }
+
+ // initialize selection log
+ if (dateAndTimes.size() <= 1) {
+ SelectInterval(it->first, it->first, itder->second, inSelection,
+ startTime, endTime);
+ if (inSelection) {
+ startTime = inputWS->run().startTime();
+ endTime = inputWS->run().endTime();
+ Kernel::SplittingInterval interval(startTime, endTime, 0);
+ splitter.push_back(interval);
+ } else {
+ throw std::runtime_error("filtered all data points. Nothing to do");
+ }
+ } else {
+ SelectInterval(it->first, next->first, itder->second, inSelection,
+ startTime, endTime);
+ }
+
+ // if its filtered using log, both iterator walk through the same values
+ // if use derivative, derivative's values are used for filtering
+ // and derivative assumed in a center of an interval
+ for (; next != dateAndTimes.end(); ++next, ++itder) {
+ if (SelectInterval(it->first, next->first, itder->second, inSelection,
+ startTime, endTime)) {
+ Kernel::SplittingInterval interval(startTime, endTime, 0);
+ splitter.push_back(interval);
+ }
+ it = next;
+ }
+ // final interval
+ if (inSelection && (endTime > startTime)) {
+ Kernel::SplittingInterval interval(startTime, endTime, 0);
+ splitter.push_back(interval);
+ }
+ } // End of USE filter log.
+
+ chop_speed = this->getAvrgLogValue(inputWS, "ChopperSpeedLog", splitter);
+ chop_speed = std::fabs(chop_speed);
+ if (chop_speed < 1.e-7) {
+ throw std::runtime_error("Chopper speed can not be zero ");
+ }
+ chop_delay =
+ std::fabs(this->getAvrgLogValue(inputWS, "ChopperDelayLog", splitter));
+
+ // process chopper delay in the units of degree (phase)
+ auto pProperty = getPLogForProperty(inputWS, "ChopperDelayLog");
+ if (!pProperty)
+ throw std::runtime_error("ChopperDelayLog has been removed from workspace "
+ "during the algorithm execution");
+ std::string units = pProperty->units();
+ // its chopper phase provided
+ if (units == "deg" ||
+ units.c_str()[0] ==
+ -80) { //<- userd in ISIS ASCII representation of o(deg)
+ chop_delay *= 1.e+6 / (360. * chop_speed); // convert in uSec
+ }
+ chop_delay += m_phase / chop_speed;
+}
+
+namespace { // namespace for lambda functions, used in validators
+
+/* former Lambda to validate if appropriate log is present in workspace
+and if it's present, it is a time-series property
+* @param prop_name -- the name of the log to check
+* @param err_presence -- core error message to return if no log found
+* @param err_type -- core error message to return if
+* log is of incorrect type
+* @param fail -- fail or warn if appropriate log is not available.
+*
+* @param result -- map to add the result of check for errors
+* if no error found the map is not modified and remains
+* empty.
+
+
+* @return -- false if all checks are fine, or true if check is
+* failed
+*/
+bool check_time_series_property(
+ const GetAllEi *algo, const API::MatrixWorkspace_sptr &inputWS,
+ const boost::shared_ptr<const Geometry::IComponent> &chopper,
+ const std::string &prop_name, const std::string &err_presence,
+ const std::string &err_type, bool fail,
+ std::map<std::string, std::string> &result) {
+
+ std::string LogName = algo->getProperty(prop_name);
+ if (boost::iequals(LogName, "Defined in IDF")) {
+ try {
+ auto theLogs = chopper->getStringParameter(prop_name);
+ if (theLogs.size() == 0) {
+ if (fail)
+ result[prop_name] = "Can not retrieve parameter " + prop_name +
+ " from the instrument definition file.";
+ return true;
+ }
+ LogName = theLogs[0];
+ } catch (...) {
+ result[prop_name] = "Can not retrieve parameter " + prop_name +
+ " from the instrument definition file.";
+ return true;
+ }
+ }
+ try {
+ Kernel::Property *pProp = inputWS->run().getProperty(LogName);
+ auto pTSProp = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(pProp);
+ if (!pTSProp) {
+ if (fail)
+ result[prop_name] = "Workspace contains " + err_type + LogName +
+ " But its type is not a timeSeries property";
+ return true;
+ }
+ } catch (std::runtime_error &) {
+ if (fail)
+ result[prop_name] = "Workspace has to contain " + err_presence + LogName;
+ return true;
+ }
+ return false;
+}
+}
+
+/**Validates if input workspace contains all necessary logs and if all
+* these logs are the logs of appropriate type
+@return list of invalid logs or empty list if no errors is found.
+*/
+std::map<std::string, std::string> GetAllEi::validateInputs() {
+
+ // Do Validation
+ std::map<std::string, std::string> result;
+
+ API::MatrixWorkspace_sptr inputWS = getProperty("Workspace");
+ if (!inputWS) {
+ result["Workspace"] = "Input workspace can not be identified";
+ return result;
+ }
+ if (!inputWS->isHistogramData()) {
+ result["Workspace"] = "Only histogram workspaces are currently supported. "
+ "Rebin input workspace first.";
+ }
+
+ specid_t specNum1 = getProperty("Monitor1SpecID");
+ try {
+ inputWS->getIndexFromSpectrumNumber(specNum1);
+ } catch (std::runtime_error &) {
+ result["Monitor1SpecID"] =
+ "Input workspace does not contain spectra with ID: " +
+ boost::lexical_cast<std::string>(specNum1);
+ }
+ specid_t specNum2 = getProperty("Monitor2SpecID");
+ try {
+ inputWS->getIndexFromSpectrumNumber(specNum2);
+ } catch (std::runtime_error &) {
+ result["Monitor2SpecID"] =
+ "Input workspace does not contain spectra with ID: " +
+ boost::lexical_cast<std::string>(specNum2);
+ }
+ // check chopper and initiate it if present (for debugging)
+ m_chopper = inputWS->getInstrument()->getComponentByName("chopper-position");
+ if (!m_chopper) {
+ result["Workspace_chopper"] =
+ " For this algorithm to work workspace has"
+ " to contain well defined 'chopper-position' component";
+ return result;
+ }
+
+ check_time_series_property(this, inputWS, m_chopper, "ChopperSpeedLog",
+ "chopper speed log with name: ",
+ "chopper speed log ", true, result);
+ check_time_series_property(
+ this, inputWS, m_chopper, "ChopperDelayLog",
+ "property related to chopper delay log with name: ", "chopper delay log ",
+ true, result);
+ bool failed = check_time_series_property(
+ this, inputWS, m_chopper, "FilterBaseLog", "filter base log named: ",
+ "filter base log: ", false, result);
+ if (failed) {
+ g_log.warning()
+ << " Can not find a log to identify good DAE operations.\n"
+ " Assuming that good operations start from experiment time=0";
+ } else {
+ this->setFilterLog(inputWS);
+ }
+ return result;
+}
+
+} // namespace Algorithms
+} // namespace Mantid
diff --git a/Framework/Algorithms/test/GetAllEiTest.h b/Framework/Algorithms/test/GetAllEiTest.h
new file mode 100644
index 0000000000000000000000000000000000000000..9dfadd5311c1761101143bc1832f02b64436d75e
--- /dev/null
+++ b/Framework/Algorithms/test/GetAllEiTest.h
@@ -0,0 +1,584 @@
+#ifndef GETALLEI_TEST_H_
+#define GETALLEI_TEST_H_
+
+#include <memory>
+#include <boost/math/special_functions/fpclassify.hpp>
+#include <cxxtest/TestSuite.h>
+#include "MantidAlgorithms/GetAllEi.h"
+#include "MantidKernel/TimeSeriesProperty.h"
+#include "MantidTestHelpers/WorkspaceCreationHelper.h"
+
+using namespace Mantid;
+using namespace Mantid::Algorithms;
+using namespace Mantid::API;
+
+class GetAllEiTester : public GetAllEi {
+public:
+ void find_chop_speed_and_delay(const API::MatrixWorkspace_sptr &inputWS,
+ double &chop_speed, double &chop_delay) {
+ GetAllEi::findChopSpeedAndDelay(inputWS, chop_speed, chop_delay);
+ }
+ void findGuessOpeningTimes(const std::pair<double, double> &TOF_range,
+ double ChopDelay, double Period,
+ std::vector<double> &guess_opening_times) {
+ GetAllEi::findGuessOpeningTimes(TOF_range, ChopDelay, Period,
+ guess_opening_times);
+ }
+ bool filterLogProvided() const { return (m_pFilterLog != NULL); }
+ double getAvrgLogValue(const API::MatrixWorkspace_sptr &inputWS,
+ const std::string &propertyName) {
+ std::vector<Kernel::SplittingInterval> splitter;
+ return GetAllEi::getAvrgLogValue(inputWS, propertyName, splitter);
+ }
+ API::MatrixWorkspace_sptr
+ buildWorkspaceToFit(const API::MatrixWorkspace_sptr &inputWS,
+ size_t &wsIndex0) {
+ return GetAllEi::buildWorkspaceToFit(inputWS, wsIndex0);
+ }
+ void findBinRanges(const MantidVec &eBins, const MantidVec &signal,
+ const std::vector<double> &guess_energies,
+ double Eresolution, std::vector<size_t> &irangeMin,
+ std::vector<size_t> &irangeMax,
+ std::vector<bool> &guessValid) {
+ GetAllEi::findBinRanges(eBins, signal, guess_energies, Eresolution,
+ irangeMin, irangeMax, guessValid);
+ }
+ void setResolution(double newResolution) {
+ this->m_max_Eresolution = newResolution;
+ }
+ size_t calcDerivativeAndCountZeros(const std::vector<double> &bins,
+ const std::vector<double> &signal,
+ std::vector<double> &deriv,
+ std::vector<double> &zeros) {
+ return GetAllEi::calcDerivativeAndCountZeros(bins, signal, deriv, zeros);
+ }
+};
+
+class GetAllEiTest : public CxxTest::TestSuite {
+public:
+ // This pair of boilerplate methods prevent the suite being created statically
+ // This means the constructor isn't called when running other tests
+ static GetAllEiTest *createSuite() { return new GetAllEiTest(); }
+ static void destroySuite(GetAllEiTest *suite) { delete suite; }
+
+ GetAllEiTest() {}
+
+public:
+ void testName() { TS_ASSERT_EQUALS(m_getAllEi.name(), "GetAllEi"); }
+
+ void testVersion() { TS_ASSERT_EQUALS(m_getAllEi.version(), 1); }
+
+ void testInit() {
+ TS_ASSERT_THROWS_NOTHING(m_getAllEi.initialize());
+ TS_ASSERT(m_getAllEi.isInitialized());
+ }
+ //
+ void test_validators_work() {
+
+ MatrixWorkspace_sptr ws = this->createTestingWS(true);
+
+ m_getAllEi.initialize();
+ m_getAllEi.setProperty("Workspace", ws);
+ m_getAllEi.setProperty("OutputWorkspace", "monitor_peaks");
+ TSM_ASSERT_THROWS(
+ "should throw runtime error on as spectra ID should be positive",
+ m_getAllEi.setProperty("Monitor1SpecID", -1), std::invalid_argument);
+
+ m_getAllEi.setProperty("Monitor1SpecID", 1);
+ m_getAllEi.setProperty("Monitor2SpecID", 2);
+ m_getAllEi.setProperty("ChopperSpeedLog", "Chopper_Speed");
+ m_getAllEi.setProperty("ChopperDelayLog", "Chopper_Delay");
+ m_getAllEi.setProperty("FilterBaseLog", "proton_charge");
+ m_getAllEi.setProperty("FilterWithDerivative", false);
+
+ TSM_ASSERT_THROWS("should throw runtime error on validation as no "
+ "appropriate logs are defined",
+ m_getAllEi.execute(), std::runtime_error);
+ auto log_messages = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("Two logs should fail", log_messages.size(), 2);
+ // add invalid property type
+ ws->mutableRun().addLogData(
+ new Kernel::PropertyWithValue<double>("Chopper_Speed", 10.));
+ auto log_messages2 = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("Two logs should fail", log_messages2.size(), 2);
+
+ TSM_ASSERT_DIFFERS("should fail for different reason ",
+ log_messages["ChopperSpeedLog"],
+ log_messages2["ChopperSpeedLog"]);
+ // add correct property type:
+ ws->mutableRun().clearLogs();
+ ws->mutableRun().addLogData(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Speed"));
+ log_messages = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("One log should fail", log_messages.size(), 1);
+ TSM_ASSERT("Filter log is not provided ", !m_getAllEi.filterLogProvided());
+ ws->mutableRun().addLogData(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Delay"));
+ ws->mutableRun().addLogData(
+ new Kernel::TimeSeriesProperty<double>("proton_charge"));
+ log_messages = m_getAllEi.validateInputs();
+
+ TSM_ASSERT_EQUALS("All logs are defined", log_messages.size(), 0);
+ TSM_ASSERT("Filter log is provided ", m_getAllEi.filterLogProvided());
+
+ m_getAllEi.setProperty("Monitor1SpecID", 3);
+ log_messages = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("Workspace should not have spectra with ID=3",
+ log_messages.size(), 1);
+ }
+ //
+ void test_get_chopper_speed() {
+
+ MatrixWorkspace_sptr ws = this->createTestingWS(true);
+
+ m_getAllEi.initialize();
+ m_getAllEi.setProperty("Workspace", ws);
+ m_getAllEi.setProperty("OutputWorkspace", "monitor_peaks");
+ m_getAllEi.setProperty("Monitor1SpecID", 1);
+ m_getAllEi.setProperty("Monitor2SpecID", 2);
+ m_getAllEi.setProperty("ChopperSpeedLog", "Chopper_Speed");
+ m_getAllEi.setProperty("ChopperDelayLog", "Chopper_Delay");
+ m_getAllEi.setProperty("FilterBaseLog", "proton_charge");
+ m_getAllEi.setProperty("FilterWithDerivative", false);
+
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> chopSpeed(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Speed"));
+ for (int i = 0; i < 10; i++) {
+ chopSpeed->addValue(Kernel::DateAndTime(10000 + 10 * i, 0), 1.);
+ }
+ for (int i = 0; i < 10; i++) {
+ chopSpeed->addValue(Kernel::DateAndTime(100 + 10 * i, 0), 10.);
+ }
+ for (int i = 0; i < 10; i++) {
+ chopSpeed->addValue(Kernel::DateAndTime(10 * i, 0), 100.);
+ }
+ ws->mutableRun().addLogData(chopSpeed.release());
+
+ // Test sort log by run time.
+ TSM_ASSERT_THROWS(
+ "Attempt to get log without start/stop time set should fail",
+ m_getAllEi.getAvrgLogValue(ws, "ChopperSpeedLog"), std::runtime_error);
+
+ ws->mutableRun().setStartAndEndTime(Kernel::DateAndTime(90, 0),
+ Kernel::DateAndTime(10000, 0));
+ double val = m_getAllEi.getAvrgLogValue(ws, "ChopperSpeedLog");
+ TS_ASSERT_DELTA(val, (10 * 10 + 100.) / 11., 1.e-6);
+
+ ws->mutableRun().setStartAndEndTime(Kernel::DateAndTime(100, 0),
+ Kernel::DateAndTime(10000, 0));
+ val = m_getAllEi.getAvrgLogValue(ws, "ChopperSpeedLog");
+ TS_ASSERT_DELTA(val, 10., 1.e-6);
+
+ // Test sort log by log.
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> chopDelay(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Delay"));
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> goodFram(
+ new Kernel::TimeSeriesProperty<double>("proton_charge"));
+
+ for (int i = 0; i < 10; i++) {
+ auto time = Kernel::DateAndTime(200 + 10 * i, 0);
+ chopDelay->addValue(time, 10.);
+ if (i < 2) {
+ goodFram->addValue(time, 1);
+ } else {
+ goodFram->addValue(time, 0);
+ }
+ }
+ for (int i = 0; i < 10; i++) {
+ auto time = Kernel::DateAndTime(100 + 10 * i, 0);
+ chopDelay->addValue(time, 0.1);
+ goodFram->addValue(time, 1);
+ }
+ for (int i = 0; i < 10; i++) {
+ auto time = Kernel::DateAndTime(10 * i, 0);
+ chopDelay->addValue(time, 1.);
+ goodFram->addValue(time, 0);
+ }
+ ws->mutableRun().addLogData(chopDelay.release());
+ ws->mutableRun().addLogData(goodFram.release());
+
+ // Run validate as this will set up property, which indicates filter log
+ // presence
+ auto errors = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("All logs are defined now", errors.size(), 0);
+
+ double chop_speed, chop_delay;
+ m_getAllEi.find_chop_speed_and_delay(ws, chop_speed, chop_delay);
+ TSM_ASSERT_DELTA("Chopper delay should have special speed ",
+ (10 * 0.1 + 20) / 12., chop_delay, 1.e-6);
+
+ goodFram.reset(new Kernel::TimeSeriesProperty<double>("proton_charge"));
+ for (int i = 0; i < 10; i++) {
+ auto time = Kernel::DateAndTime(100 + 10 * i, 0);
+ goodFram->addValue(time, 1);
+ }
+
+ ws->mutableRun().addProperty(goodFram.release(), true);
+ errors = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("All logs are defined now", errors.size(), 0);
+
+ m_getAllEi.find_chop_speed_and_delay(ws, chop_speed, chop_delay);
+ TSM_ASSERT_DELTA("Chopper delay should have special speed", 0.1, chop_delay,
+ 1.e-6);
+ }
+ void test_get_chopper_speed_filter_derivative() {
+
+ MatrixWorkspace_sptr ws = this->createTestingWS(true);
+
+ m_getAllEi.initialize();
+ m_getAllEi.setProperty("Workspace", ws);
+ m_getAllEi.setProperty("OutputWorkspace", "monitor_peaks");
+ m_getAllEi.setProperty("Monitor1SpecID", 1);
+ m_getAllEi.setProperty("Monitor2SpecID", 2);
+ m_getAllEi.setProperty("ChopperSpeedLog", "Chopper_Speed");
+ m_getAllEi.setProperty("ChopperDelayLog", "Chopper_Delay");
+ m_getAllEi.setProperty("FilterBaseLog", "proton_charge");
+ m_getAllEi.setProperty("FilterWithDerivative", true);
+
+ // Test select log by log derivative
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> chopDelay(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Delay"));
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> chopSpeed(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Speed"));
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> protCharge(
+ new Kernel::TimeSeriesProperty<double>("proton_charge"));
+
+ double gf(0);
+ for (int i = 0; i < 50; i++) {
+ auto time = Kernel::DateAndTime(10 * i, 0);
+ if (i > 10 && i < 20) {
+ chopDelay->addValue(time, 100.);
+ chopSpeed->addValue(time, 0.);
+ protCharge->addValue(time, gf);
+ } else {
+ chopDelay->addValue(time, 10.);
+ chopSpeed->addValue(time, 50.);
+ protCharge->addValue(time, gf);
+ gf++;
+ }
+ }
+ ws->mutableRun().addLogData(chopSpeed.release());
+ ws->mutableRun().addLogData(chopDelay.release());
+ ws->mutableRun().addLogData(protCharge.release());
+
+ // Run validate as this will set up property, which indicates filter log
+ // presence
+ auto errors = m_getAllEi.validateInputs();
+ TSM_ASSERT_EQUALS("All logs are defined now", errors.size(), 0);
+
+ double chop_speed, chop_delay;
+ m_getAllEi.find_chop_speed_and_delay(ws, chop_speed, chop_delay);
+ TSM_ASSERT_DELTA("Chopper delay should have defined value ", 10.,
+ chop_delay, 1.e-6);
+ TSM_ASSERT_DELTA("Chopper speed should have defined speed", 50., chop_speed,
+ 1.e-6);
+ }
+
+ void test_guess_opening_times() {
+
+ std::pair<double, double> TOF_range(5, 100);
+ double t0(6), Period(10);
+ std::vector<double> guess_tof;
+ m_getAllEi.findGuessOpeningTimes(TOF_range, t0, Period, guess_tof);
+ TSM_ASSERT_EQUALS("should have 10 periods within the specified interval",
+ guess_tof.size(), 10);
+
+ guess_tof.resize(0);
+ t0 = TOF_range.first;
+ m_getAllEi.findGuessOpeningTimes(TOF_range, t0, Period, guess_tof);
+ TSM_ASSERT_EQUALS(
+ "Still should be 10 periods within the specified interval",
+ guess_tof.size(), 10);
+
+ t0 = TOF_range.second;
+ TSM_ASSERT_THROWS(
+ "Should throw out of range",
+ m_getAllEi.findGuessOpeningTimes(TOF_range, t0, Period, guess_tof),
+ std::runtime_error);
+
+ t0 = 1;
+ guess_tof.resize(0);
+ m_getAllEi.findGuessOpeningTimes(TOF_range, t0, Period, guess_tof);
+ TSM_ASSERT_EQUALS(" should be 9 periods within the specified interval",
+ guess_tof.size(), 9);
+
+ guess_tof.resize(0);
+ t0 = 21;
+ TOF_range.first = 20;
+ m_getAllEi.findGuessOpeningTimes(TOF_range, t0, Period, guess_tof);
+ TSM_ASSERT_EQUALS(" should be 8 periods within the specified interval",
+ guess_tof.size(), 8);
+ }
+ //
+ void test_internalWS_to_fit() {
+ Mantid::DataObjects::Workspace2D_sptr tws =
+ WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(5, 100,
+ true);
+ auto det1 = tws->getDetector(0);
+ auto det2 = tws->getDetector(4);
+ auto spec1 = tws->getSpectrum(0);
+ auto spec2 = tws->getSpectrum(4);
+ auto detID1 = spec1->getDetectorIDs();
+ auto detID2 = spec2->getDetectorIDs();
+
+ m_getAllEi.initialize();
+ m_getAllEi.setProperty("Workspace", tws);
+ m_getAllEi.setProperty("OutputWorkspace", "monitor_peaks");
+ m_getAllEi.setProperty("Monitor1SpecID", 1);
+ m_getAllEi.setProperty("Monitor2SpecID", 5);
+
+ size_t wsIndex0;
+ auto wws = m_getAllEi.buildWorkspaceToFit(tws, wsIndex0);
+
+ auto det1p = wws->getDetector(0);
+ auto det2p = wws->getDetector(1);
+ TSM_ASSERT_EQUALS("should be the same first detector position",
+ det1p->getRelativePos(), det1->getRelativePos());
+ TSM_ASSERT_EQUALS("should be the same second detector position",
+ det2p->getRelativePos(), det2->getRelativePos());
+
+ TSM_ASSERT_EQUALS("Detector's ID for the first spectrum and new workspace "
+ "should coincide",
+ *(detID1.begin()),
+ (*wws->getSpectrum(0)->getDetectorIDs().begin()));
+ TSM_ASSERT_EQUALS("Detector's ID for the second spectrum and new workspace "
+ "should coincide",
+ *(detID2.begin()),
+ (*wws->getSpectrum(1)->getDetectorIDs().begin()));
+ auto pSpec1 = wws->getSpectrum(0);
+ auto pSpec2 = wws->getSpectrum(1);
+ auto Xsp1 = pSpec1->dataX();
+ auto Xsp2 = pSpec2->dataX();
+ size_t nSpectra = Xsp2.size();
+ TS_ASSERT_EQUALS(nSpectra, 101);
+ TS_ASSERT(boost::math::isinf(Xsp1[nSpectra - 1]));
+ TS_ASSERT(boost::math::isinf(Xsp2[nSpectra - 1]));
+
+ // for(size_t i=0;i<Xsp1.size();i++){
+ // TS_ASSERT_DELTA(Xsp1[i],Xsp2[i],1.e-6);
+ //}
+ }
+ void test_calcDerivative() {
+ double sig[] = {1, 2, 3, 4, 5, 6};
+ std::vector<double> signal(sig, sig + sizeof(sig) / sizeof(double));
+ double bin[] = {2, 3, 4, 5, 6, 7, 8};
+ std::vector<double> bins(bin, bin + sizeof(bin) / sizeof(double));
+ std::vector<double> zeros;
+
+ std::vector<double> deriv;
+ size_t nZer =
+ m_getAllEi.calcDerivativeAndCountZeros(bins, signal, deriv, zeros);
+ TS_ASSERT_EQUALS(nZer, 0);
+ TS_ASSERT_DELTA(deriv[0], deriv[1], 1.e-9);
+ TS_ASSERT_DELTA(deriv[0], deriv[5], 1.e-9);
+ TS_ASSERT_DELTA(deriv[0], deriv[2], 1.e-9);
+ TS_ASSERT_DELTA(deriv[0], 1., 1.e-9);
+
+ double bin1[] = {0, 1, 3, 6, 10, 15, 21};
+ std::vector<double> bins1(bin1, bin1 + sizeof(bin1) / sizeof(double));
+ nZer = m_getAllEi.calcDerivativeAndCountZeros(bins1, signal, deriv, zeros);
+ TS_ASSERT_EQUALS(nZer, 0);
+ TS_ASSERT_DELTA(deriv[0], deriv[1], 1.e-9);
+ TS_ASSERT_DELTA(deriv[0], deriv[5], 1.e-9);
+ TS_ASSERT_DELTA(deriv[0], deriv[2], 1.e-9);
+ TS_ASSERT_DELTA(deriv[0], 0, 1.e-9);
+
+ bins.resize(101);
+ signal.resize(100);
+ for (size_t i = 0; i < 101; i++) {
+ bins[i] = double(i) * 0.1;
+ }
+ for (size_t i = 0; i < 100; i++) {
+ signal[i] = std::sin(0.5 * (bins[i] + bins[i + 1]));
+ }
+ nZer = m_getAllEi.calcDerivativeAndCountZeros(bins, signal, deriv, zeros);
+ TS_ASSERT_EQUALS(nZer, 3);
+ for (size_t i = 0; i < 99; i++) { // intentionally left boundary point --
+ // its accuracy is much lower
+ TSM_ASSERT_DELTA("At i=" + boost::lexical_cast<std::string>(i), deriv[i],
+ 10. * std::cos(0.5 * (bins[i] + bins[i + 1])), 1.e-1);
+ }
+ TS_ASSERT_DELTA(zeros[0], 1.55, 1.e-3);
+ TS_ASSERT_DELTA(zeros[1], 4.65, 1.e-3);
+ TS_ASSERT_DELTA(zeros[2], 7.85, 1.e-3);
+ }
+ void test_binRanges() {
+ std::vector<size_t> bin_min, bin_max, zeros;
+ // Index 0 1 2 3 4 5 6 7 8 9 10 11 12 13
+ double debin[] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15};
+ std::vector<double> ebin(debin, debin + sizeof(debin) / sizeof(double));
+ // Not yet supported in VC 2012
+ // std::vector<double> ebin={1,2,3,4,5,6,7,8,9,10,11,12,13,15};
+ // Ind 0 1 2 3 4 5 6 7 8 9 10 11 12 13
+ double sig[] = {0, 0, 0, 3, 0, 0, 4, 0, 0, 0, 11, 0, 0};
+ std::vector<double> signal(sig, sig + sizeof(sig) / sizeof(double));
+
+ double dguess[] = {1, 6, 10, 12};
+ std::vector<double> guess(dguess, dguess + sizeof(dguess) / sizeof(double));
+ std::vector<bool> guessValid;
+
+ m_getAllEi.findBinRanges(ebin, signal, guess, 0.1, bin_min, bin_max,
+ guessValid);
+
+ TS_ASSERT_EQUALS(bin_min.size(), 2)
+ TS_ASSERT_EQUALS(bin_max.size(), 2)
+ TS_ASSERT_EQUALS(guessValid.size(), 4)
+ TS_ASSERT_EQUALS(bin_min[0], 4)
+ TS_ASSERT_EQUALS(bin_max[0], 9)
+ TS_ASSERT_EQUALS(bin_min[1], 7)
+ TS_ASSERT_EQUALS(bin_max[1], 13)
+
+ signal[10] = 0;
+ signal[11] = 11;
+ guess[1] = 3;
+ guess[2] = 6;
+ guess[3] = 11;
+ m_getAllEi.findBinRanges(ebin, signal, guess, 0.01, bin_min, bin_max,
+ guessValid);
+ TS_ASSERT_EQUALS(bin_min.size(), 3)
+ TS_ASSERT_EQUALS(bin_max.size(), 3)
+ TS_ASSERT_EQUALS(guessValid.size(), 4)
+
+ TS_ASSERT_EQUALS(bin_min[0], 3);
+ TS_ASSERT_EQUALS(bin_max[0], 4);
+ TS_ASSERT(guessValid[1]);
+
+ TS_ASSERT_EQUALS(bin_min[1], 6);
+ TS_ASSERT_EQUALS(bin_max[1], 7);
+ TS_ASSERT(guessValid[2]);
+
+ TS_ASSERT_EQUALS(bin_min[2], 11);
+ TS_ASSERT_EQUALS(bin_max[2], 12);
+ TS_ASSERT(guessValid[3]);
+
+ TS_ASSERT(!guessValid[0]);
+ }
+
+ void test_getAllEi() {
+ auto ws = createTestingWS();
+
+ m_getAllEi.initialize();
+ m_getAllEi.setProperty("Workspace", ws);
+ m_getAllEi.setProperty("OutputWorkspace", "monitor_peaks");
+ m_getAllEi.setProperty("Monitor1SpecID", 1);
+ m_getAllEi.setProperty("Monitor2SpecID", 2);
+ m_getAllEi.setProperty("ChopperSpeedLog", "Chopper_Speed");
+ m_getAllEi.setProperty("ChopperDelayLog", "Chopper_Delay");
+ m_getAllEi.setProperty("FilterBaseLog", "is_running");
+ m_getAllEi.setProperty("FilterWithDerivative", false);
+ m_getAllEi.setProperty("OutputWorkspace", "allEiWs");
+
+ TS_ASSERT_THROWS_NOTHING(m_getAllEi.execute());
+ API::MatrixWorkspace_sptr out_ws;
+ TS_ASSERT_THROWS_NOTHING(
+ out_ws = API::AnalysisDataService::Instance()
+ .retrieveWS<API::MatrixWorkspace>("allEiWs"));
+
+ TSM_ASSERT("Should be able to retrieve workspace", out_ws);
+ auto wso = dynamic_cast<DataObjects::Workspace2D *>(out_ws.get());
+ TS_ASSERT(wso);
+ if (!wso)
+ return;
+
+ auto &x = wso->dataX(0);
+ TSM_ASSERT_EQUALS("Second peak should be filtered by monitor ranges",
+ x.size(), 1);
+ TS_ASSERT_DELTA(x[0], 134.316, 1.e-3)
+ }
+
+private:
+ GetAllEiTester m_getAllEi;
+
+ DataObjects::Workspace2D_sptr createTestingWS(bool noLogs = false) {
+ double delay(2000), chopSpeed(100), inital_chop_phase(-3000);
+ auto ws = WorkspaceCreationHelper::create2DWorkspaceWithFullInstrument(
+ 2, 1000, true);
+ auto pInstrument = ws->getInstrument();
+ auto chopper = pInstrument->getComponentByName("chopper-position");
+
+ // add chopper parameters
+ auto ¶mMap = ws->instrumentParameters();
+ const std::string description(
+ "The initial rotation phase of the disk used to calculate the time"
+ " for neutrons arriving at the chopper according to the formula time = "
+ "delay + initial_phase/Speed");
+ paramMap.add<double>("double", chopper.get(), "initial_phase",
+ inital_chop_phase, &description);
+ paramMap.add<std::string>("string", chopper.get(), "ChopperDelayLog",
+ "fermi_delay");
+ paramMap.add<std::string>("string", chopper.get(), "ChopperSpeedLog",
+ "fermi_speed");
+ paramMap.add<std::string>("string", chopper.get(), "FilterBaseLog",
+ "is_running");
+ paramMap.add<bool>("bool", chopper.get(), "filter_with_derivative", false);
+
+ // test instrument parameters (obtained from workspace):
+ auto moderator = pInstrument->getSource();
+ auto detector1 = ws->getDetector(0);
+ auto detector2 = ws->getDetector(1);
+ double l_chop = chopper->getDistance(*moderator);
+ double l_mon1 = detector1->getDistance(*moderator);
+ double l_mon2 = detector2->getDistance(*moderator);
+ //,l_mon1(20-9),l_mon2(20-2);
+ double t_chop(delay + inital_chop_phase / chopSpeed);
+ double Period =
+ (0.5 * 1.e+6) / chopSpeed; // 0.5 because some choppers open twice.
+ auto &x = ws->dataX(0);
+ for (size_t i = 0; i < x.size(); i++) {
+ x[i] = 5 + double(i) * 10;
+ }
+ // signal at first monitor
+ double t1 = t_chop * l_mon1 / l_chop;
+ double t2 = (t_chop + Period) * l_mon1 / l_chop;
+ {
+ auto &y = ws->dataY(0);
+ for (size_t i = 0; i < y.size(); i++) {
+ double t = 0.5 * (x[i] + x[i + 1]);
+ double tm1 = t - t1;
+ double tm2 = t - t2;
+ y[i] = (10000 * std::exp(-tm1 * tm1 / 1000.) +
+ 20000 * std::exp(-tm2 * tm2 / 1000.));
+ // std::cout<<"t="<<t<<" signal="<<y[i]<<" ind="<<i<<std::endl;
+ }
+ }
+ // signal at second monitor
+ t1 = t_chop * l_mon2 / l_chop;
+ t2 = (t_chop + Period) * l_mon2 / l_chop;
+ {
+ auto &y = ws->dataY(1);
+ for (size_t i = 0; i < y.size(); i++) {
+ double t = 0.5 * (x[i] + x[i + 1]);
+ double tm1 = t - t1;
+ double tm2 = t - t2;
+ y[i] = (100 * std::exp(-tm1 * tm1 / 1000.) +
+ 200 * std::exp(-tm2 * tm2 / 1000.));
+ // std::cout<<"t="<<t<<" signal="<<y[i]<<" ind="<<i<<std::endl;
+ }
+ }
+
+ if (noLogs)
+ return ws;
+
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> chopDelayLog(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Delay"));
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> chopSpeedLog(
+ new Kernel::TimeSeriesProperty<double>("Chopper_Speed"));
+ std::unique_ptr<Kernel::TimeSeriesProperty<double>> isRunning(
+ new Kernel::TimeSeriesProperty<double>("is_running"));
+
+ for (int i = 0; i < 10; i++) {
+ auto time = Kernel::DateAndTime(10 * i, 0);
+ chopDelayLog->addValue(time, delay);
+ chopSpeedLog->addValue(time, chopSpeed);
+ isRunning->addValue(time, 1.);
+ }
+
+ ws->mutableRun().addLogData(chopSpeedLog.release());
+ ws->mutableRun().addLogData(chopDelayLog.release());
+ ws->mutableRun().addLogData(isRunning.release());
+
+ return ws;
+ }
+};
+
+#endif
diff --git a/Framework/CurveFitting/src/SimplexMinimizer.cpp b/Framework/CurveFitting/src/SimplexMinimizer.cpp
index b0c5ba671041eaada2fdb2979a8abb0249bd05dc..694b344296ef6cb9944ff1f1932b6950cc8bd19e 100644
--- a/Framework/CurveFitting/src/SimplexMinimizer.cpp
+++ b/Framework/CurveFitting/src/SimplexMinimizer.cpp
@@ -55,6 +55,9 @@ void SimplexMinimizer::initialize(API::ICostFunction_sptr function, size_t) {
size_t np = function->nParams();
// step size for simplex
m_simplexStepSize = gsl_vector_alloc(np);
+ if (m_simplexStepSize == NULL) {
+ throw std::runtime_error("Simplex step size initialized to NULL pointer");
+ }
gsl_vector_set_all(m_simplexStepSize, m_size);
// setup simplex container
diff --git a/Framework/Geometry/src/Instrument.cpp b/Framework/Geometry/src/Instrument.cpp
index 8841dfdceac65fa7349966abb79781c368247339..0990d31e19b2e3d79b816f62485c16a9f91b945f 100644
--- a/Framework/Geometry/src/Instrument.cpp
+++ b/Framework/Geometry/src/Instrument.cpp
@@ -1258,6 +1258,11 @@ Instrument::ContainsState Instrument::containsRectDetectors() const {
if (detector && detector->isMonitor())
continue;
+ // skip choppers HACK!
+ if (comp->getName() == "chopper-position") {
+ continue;
+ }
+
if (dynamic_cast<const RectangularDetector *>(comp.get())) {
if (!foundRect)
foundRect = true;
diff --git a/Framework/Kernel/inc/MantidKernel/TimeSeriesProperty.h b/Framework/Kernel/inc/MantidKernel/TimeSeriesProperty.h
index 46177e23dcdaaccdd8236528af0aacc35c87d13a..cf1df3f3df0a80b93e127d1c0f38b4778ef017c9 100644
--- a/Framework/Kernel/inc/MantidKernel/TimeSeriesProperty.h
+++ b/Framework/Kernel/inc/MantidKernel/TimeSeriesProperty.h
@@ -116,6 +116,11 @@ public:
virtual ~TimeSeriesProperty();
/// "Virtual" copy constructor
TimeSeriesProperty<TYPE> *clone() const;
+ //
+ /// Return time series property, containing time derivative of current
+ /// property
+ std::unique_ptr<TimeSeriesProperty<double>> getDerivative() const;
+
/// "Virtual" copy constructor with a time shift in seconds
virtual Property *cloneWithTimeShift(const double timeShift) const;
/// Return the memory used by the property, in bytes
@@ -273,6 +278,11 @@ public:
/// Stringize the property
std::string toString() const;
+ /**Reserve memory for efficient adding values to existing property
+ * makes sense only when you have reasonably precise estimate of the
+ * total size you'll need easily available in advance. */
+ void reserve(size_t size) { m_values.reserve(size); };
+
private:
/// Sort the property into increasing times
void sort() const;
diff --git a/Framework/Kernel/inc/MantidKernel/VectorHelper.h b/Framework/Kernel/inc/MantidKernel/VectorHelper.h
index b1643b23acca49c8d9be328558d7e745338d3f79..489b6d1864b0407af065d5b6301f97c7d0980607 100644
--- a/Framework/Kernel/inc/MantidKernel/VectorHelper.h
+++ b/Framework/Kernel/inc/MantidKernel/VectorHelper.h
@@ -59,7 +59,7 @@ rebinHistogram(const std::vector<double> &xold, const std::vector<double> &yold,
std::vector<double> &ynew, std::vector<double> &enew,
bool addition);
-/// Convert an array of bin boundaries to bin centre values.
+/// Convert an array of bin boundaries to bin center values.
void MANTID_KERNEL_DLL convertToBinCentre(const std::vector<double> &bin_edges,
std::vector<double> &bin_centres);
@@ -81,6 +81,15 @@ MANTID_KERNEL_DLL int getBinIndex(const std::vector<double> &bins,
MANTID_KERNEL_DLL void linearlyInterpolateY(const std::vector<double> &x,
std::vector<double> &y,
const double stepSize);
+// Do running average of input vector within specified range, considering
+// heterogeneous bin-boundaries
+// if such boundaries are provided
+MANTID_KERNEL_DLL void
+smoothInRange(const std::vector<double> &input, std::vector<double> &output,
+ double avrgInterval,
+ std::vector<double> const *const binBoundaris = NULL,
+ size_t startIndex = 0, size_t endIndex = 0,
+ std::vector<double> *const outputBinBoundaries = NULL);
//-------------------------------------------------------------------------------------
/** Return the length of the vector (in the physical sense),
@@ -102,7 +111,7 @@ template <typename T>
T scalar_prod(const std::vector<T> &v1, const std::vector<T> &v2) {
if (v1.size() != v2.size())
throw std::invalid_argument(" scalar product is defined only for the "
- "vectors of the equivalient length");
+ "vectors of the equivalent length");
T total = 0;
for (size_t i = 0; i < v1.size(); i++)
total += v1[i] * v2[i];
@@ -153,7 +162,7 @@ template <class T> struct SumGaussError : public std::binary_function<T, T, T> {
};
/**
- * Functor to deal with the increase in the error when adding (or substracting)
+ * Functor to deal with the increase in the error when adding (or subtracting)
* a number of counts.
* More generally add errors in quadrature using the square of one of the errors
* (variance = error^2)
diff --git a/Framework/Kernel/src/TimeSeriesProperty.cpp b/Framework/Kernel/src/TimeSeriesProperty.cpp
index 3cc831022f25f516b30603ebccbaca119f83f68f..f1a10e4b434adeb6455e9f99648ed6ac8673f76f 100644
--- a/Framework/Kernel/src/TimeSeriesProperty.cpp
+++ b/Framework/Kernel/src/TimeSeriesProperty.cpp
@@ -57,6 +57,54 @@ TimeSeriesProperty<TYPE>::cloneWithTimeShift(const double timeShift) const {
return timeSeriesProperty;
}
+/** Return time series property, containing time derivative of current property.
+* The property itself and the returned time derivative become sorted by time and
+* the derivative is calculated in seconds^-1.
+* (e.g. dValue/dT where dT=t2-t1 is time difference in seconds
+* for subsequent time readings and dValue=Val1-Val2 is difference in
+* subsequent values)
+*
+*/
+template <typename TYPE>
+std::unique_ptr<TimeSeriesProperty<double>>
+TimeSeriesProperty<TYPE>::getDerivative() const {
+
+ if (this->m_values.size() < 2) {
+ throw std::runtime_error("Derivative is not defined for a time-series "
+ "property with less then two values");
+ }
+
+ this->sort();
+ auto it = this->m_values.begin();
+ int64_t t0 = it->time().totalNanoseconds();
+ TYPE v0 = it->value();
+
+ it++;
+ auto timeSeriesDeriv = std::unique_ptr<TimeSeriesProperty<double>>(
+ new TimeSeriesProperty<double>(this->name() + "_derivative"));
+ timeSeriesDeriv->reserve(this->m_values.size() - 1);
+ for (; it != m_values.end(); it++) {
+ TYPE v1 = it->value();
+ int64_t t1 = it->time().totalNanoseconds();
+ if (t1 != t0) {
+ double deriv = 1.e+9 * (double(v1 - v0) / double(t1 - t0));
+ int64_t tm = static_cast<int64_t>((t1 + t0) / 2);
+ timeSeriesDeriv->addValue(Kernel::DateAndTime(tm), deriv);
+ }
+ t0 = t1;
+ v0 = v1;
+ }
+ return timeSeriesDeriv;
+}
+/** time series derivative specialization for string type */
+template <>
+std::unique_ptr<TimeSeriesProperty<double>>
+TimeSeriesProperty<std::string>::getDerivative() const {
+ throw std::runtime_error(
+ "Time series property derivative is not defined for strings");
+ // return nullptr;
+}
+
/**
* Return the memory used by the property, in bytes
* */
diff --git a/Framework/Kernel/src/VectorHelper.cpp b/Framework/Kernel/src/VectorHelper.cpp
index 61aedadfc83fefb11da536e3229170dc378d050e..52144d225aa4f454563c018cefb3c214ad817445 100644
--- a/Framework/Kernel/src/VectorHelper.cpp
+++ b/Framework/Kernel/src/VectorHelper.cpp
@@ -533,6 +533,164 @@ void linearlyInterpolateY(const std::vector<double> &x, std::vector<double> &y,
step++;
}
}
+namespace {
+/** internal function converted from Lambda to identify interval around
+* specified point and run average around this point
+*
+*@param index -- index to average around
+*@param startIndex -- index in the array of data (input to start average
+* from) should be: index>=startIndex>=0
+*@param endIndex -- index in the array of data (input to end average at)
+* should be: index<=endIndex<=input.size()
+*@param halfWidth -- half width of the interval to integrate.
+*@param input -- vector of input signal
+*@param binBndrs -- pointer to vector of bin boundaries or NULL pointer.
+ */
+double runAverage(size_t index, size_t startIndex, size_t endIndex,
+ const double halfWidth, const std::vector<double> &input,
+ std::vector<double> const *const binBndrs) {
+
+ size_t iStart, iEnd;
+ double weight0(0), weight1(0), start, end;
+ //
+ if (binBndrs) {
+ // identify initial and final bins to
+ // integrate over. Notice the difference
+ // between start and end bin and shift of
+ // the interpolating function into the center
+ // of each bin
+ auto &rBndrs = *binBndrs;
+ // bin0 = binBndrs->operator[](index + 1) - binBndrs->operator[](index);
+
+ double binC = 0.5 * (rBndrs[index + 1] + rBndrs[index]);
+ start = binC - halfWidth;
+ end = binC + halfWidth;
+ if (start <= rBndrs[startIndex]) {
+ iStart = startIndex;
+ start = rBndrs[iStart];
+ } else {
+ iStart = getBinIndex(*binBndrs, start);
+ weight0 =
+ (rBndrs[iStart + 1] - start) / (rBndrs[iStart + 1] - rBndrs[iStart]);
+ iStart++;
+ }
+ if (end >= rBndrs[endIndex]) {
+ iEnd = endIndex; // the signal defined up to i<iEnd
+ end = rBndrs[endIndex];
+ } else {
+ iEnd = getBinIndex(*binBndrs, end);
+ weight1 = (end - rBndrs[iEnd]) / (rBndrs[iEnd + 1] - rBndrs[iEnd]);
+ }
+ if (iStart > iEnd) { // start and end get into the same bin
+ weight1 = 0;
+ weight0 = (end - start) / (rBndrs[iStart] - rBndrs[iStart - 1]);
+ }
+ } else { // integer indexes and functions defined in the bin centers
+ auto iHalfWidth = static_cast<size_t>(halfWidth);
+ iStart = index - iHalfWidth;
+ if (startIndex + iHalfWidth > index)
+ iStart = startIndex;
+ iEnd = index + iHalfWidth;
+ if (iEnd > endIndex)
+ iEnd = endIndex;
+ }
+
+ double avrg = 0;
+ size_t ic = 0;
+ for (size_t j = iStart; j < iEnd; j++) {
+ avrg += input[j];
+ ic++;
+ }
+ if (binBndrs) { // add values at edges
+ if (iStart != startIndex)
+ avrg += input[iStart - 1] * weight0;
+ if (iEnd != endIndex)
+ avrg += input[iEnd] * weight1;
+
+ return avrg / (end - start);
+ } else {
+ return avrg / double(ic);
+ }
+}
+}
+/** Basic running average of input vector within specified range, considering
+* variable bin-boundaries if such boundaries are provided.
+* The algorithm performs trapezium integration, so some peak shift
+* related to the first derivative of the integrated function can be observed.
+*
+* @param input:: input vector to smooth
+* @param output:: resulting vector (can not coincide with input)
+* @param avrgInterval:: the interval to average function in.
+* the function is averaged within +-0.5*avrgInterval
+* @param binBndrs :: pointer to the vector, containing bin boundaries.
+* If provided, its length has to be input.size()+1,
+* if not, equal size bins of size 1 are assumed,
+* so avrgInterval becomes the number of points
+* to average over. Bin boundaries array have to
+* increase and can not contain equal boundaries.
+* @param startIndex:: if provided, its start index to run averaging from.
+* if not, averaging starts from the index 0
+* @param endIndex :: final index to run average to, if provided. If
+* not, or higher then number of elements in input array,
+* averaging is performed to the end point of the input
+* array
+* @param outBins :: if present, pointer to a vector to return
+* bin boundaries for output array.
+*/
+void smoothInRange(const std::vector<double> &input,
+ std::vector<double> &output, const double avrgInterval,
+ std::vector<double> const *const binBndrs, size_t startIndex,
+ size_t endIndex, std::vector<double> *const outBins) {
+
+ if (endIndex == 0)
+ endIndex = input.size();
+ if (endIndex > input.size())
+ endIndex = input.size();
+
+ if (endIndex <= startIndex) {
+ output.resize(0);
+ return;
+ }
+
+ size_t max_size = input.size();
+ if (binBndrs) {
+ if (binBndrs->size() != max_size + 1) {
+ throw std::invalid_argument(
+ "Array of bin boundaries, "
+ "if present, have to be one bigger then the input array");
+ }
+ }
+
+ size_t length = endIndex - startIndex;
+ output.resize(length);
+
+ double halfWidth = avrgInterval / 2;
+ if (!binBndrs) {
+ if (std::floor(halfWidth) * 2 - avrgInterval > 1.e-6) {
+ halfWidth = std::floor(halfWidth) + 1;
+ }
+ }
+
+ if (outBins)
+ outBins->resize(length + 1);
+
+ // Run averaging
+ double binSize = 1;
+ for (size_t i = startIndex; i < endIndex; i++) {
+ if (binBndrs) {
+ binSize = binBndrs->operator[](i + 1) - binBndrs->operator[](i);
+ }
+ output[i - startIndex] =
+ runAverage(i, startIndex, endIndex, halfWidth, input, binBndrs) *
+ binSize;
+ if (outBins) {
+ outBins->operator[](i - startIndex) = binBndrs->operator[](i);
+ }
+ }
+ if (outBins) {
+ outBins->operator[](endIndex - startIndex) = binBndrs->operator[](endIndex);
+ }
+}
/// Declare all version of this
template DLLExport std::vector<int32_t>
diff --git a/Framework/Kernel/test/TimeSeriesPropertyTest.h b/Framework/Kernel/test/TimeSeriesPropertyTest.h
index 284bfb9bc6f55d75e19630a1788b1839afbb4727..87a1480b405e3035e1babbded9d0a4469e34a2c5 100644
--- a/Framework/Kernel/test/TimeSeriesPropertyTest.h
+++ b/Framework/Kernel/test/TimeSeriesPropertyTest.h
@@ -122,7 +122,36 @@ public:
TS_ASSERT_EQUALS(twoVals[1], threeVals[2]);
TS_ASSERT_EQUALS(newVal, threeVals[1]);
}
+ void test_GetDerivative() {
+ dProp->addValue("2007-11-30T16:17:10", 10);
+ dProp->addValue("2007-11-30T16:17:12", 12);
+ dProp->addValue("2007-11-30T16:17:01", 01);
+ dProp->addValue("2007-11-30T16:17:05", 05);
+ auto derProp = dProp->getDerivative();
+ TS_ASSERT(dynamic_cast<TimeSeriesProperty<double> *>(derProp.get()))
+
+ TS_ASSERT_EQUALS(derProp->size(), 3);
+ auto derValues = derProp->valuesAsVector();
+
+ TS_ASSERT_EQUALS(derValues[0], 1);
+ TS_ASSERT_EQUALS(derValues[1], 1);
+ TS_ASSERT_EQUALS(derValues[2], 1);
+
+ TSM_ASSERT_THROWS("derivative undefined for string property",
+ sProp->getDerivative(), std::runtime_error);
+
+ iProp->addValue("2007-11-30T16:17:10", 10);
+ TSM_ASSERT_THROWS(
+ "derivative undefined for property with less then 2 values",
+ iProp->getDerivative(), std::runtime_error);
+ iProp->addValue("2007-11-30T16:17:12", 12);
+
+ derProp = iProp->getDerivative();
+ TS_ASSERT_EQUALS(derProp->size(), 1);
+ derValues = derProp->valuesAsVector();
+ TS_ASSERT_EQUALS(derValues[0], 1);
+ }
void test_timesAsVector() {
TimeSeriesProperty<double> *p =
new TimeSeriesProperty<double>("doubleProp");
diff --git a/Framework/Kernel/test/VectorHelperTest.h b/Framework/Kernel/test/VectorHelperTest.h
index 16f9f4b86a2f79e2fe41a07dc55323e44066e000..36ac636a8f3644be183ca269644490a15f7f741e 100644
--- a/Framework/Kernel/test/VectorHelperTest.h
+++ b/Framework/Kernel/test/VectorHelperTest.h
@@ -5,6 +5,7 @@
#include "MantidKernel/Timer.h"
#include "MantidKernel/VectorHelper.h"
#include <cxxtest/TestSuite.h>
+#include <cstdlib>
#include <boost/assign/list_of.hpp>
using namespace Mantid::Kernel;
@@ -266,6 +267,121 @@ public:
index = VectorHelper::getBinIndex(m_test_bins, testValue));
TS_ASSERT_EQUALS(index, 2);
}
+ void test_RunningAveraging() {
+ double id[] = {1, 2, 3, 4, 5, 6};
+ std::vector<double> inputData(id, id + sizeof(id) / sizeof(double));
+ double ib[] = {0, 1, 2, 3, 4, 5};
+ std::vector<double> inputBoundaries(ib, ib + sizeof(ib) / sizeof(double));
+
+ std::vector<double> output;
+ TS_ASSERT_THROWS(
+ VectorHelper::smoothInRange(inputData, output, 6, &inputBoundaries),
+ std::invalid_argument);
+ inputBoundaries.push_back(6);
+ VectorHelper::smoothInRange(inputData, output, 6, &inputBoundaries);
+
+ TS_ASSERT_DELTA(output[1] - output[0], 0.492, 1.e-3);
+ TS_ASSERT_DELTA(output[3] - output[2], 0.4545, 1.e-3);
+ TS_ASSERT_DELTA(output[5] - output[4], 0.492, 1.e-3);
+ inputBoundaries[1] = 1;
+ inputBoundaries[2] = 3;
+ inputBoundaries[3] = 6;
+ inputBoundaries[4] = 10;
+ inputBoundaries[5] = 15;
+ inputBoundaries[6] = 21;
+ VectorHelper::smoothInRange(inputData, output, 6, &inputBoundaries);
+ TS_ASSERT_DELTA(output[2], 3, 1.e-8);
+ TS_ASSERT_DELTA(output[0], 1, 1.e-8);
+ TS_ASSERT_DELTA(output[5], 6, 1.e-8);
+
+ std::vector<double> out_bins;
+ VectorHelper::smoothInRange(inputData, output, 3, &inputBoundaries, 1, 5,
+ &out_bins);
+ TS_ASSERT_EQUALS(output.size(), 4);
+ TS_ASSERT_DELTA(output[1], 3, 1.e-8);
+ }
+
+ void test_Smooth_keeps_peakPosition() {
+
+ std::vector<double> output;
+ std::vector<double> inputBoundaries(21);
+ inputBoundaries[0] = 0;
+ double step(1);
+ for (size_t i = 1; i < 21; i++) {
+ inputBoundaries[i] = inputBoundaries[i - 1] + step;
+ step *= 1.1;
+ }
+ double norm = 100 / inputBoundaries[20];
+ for (size_t i = 0; i < 21; i++) {
+ inputBoundaries[i] *= norm;
+ }
+
+ std::vector<double> inputData(20);
+ for (size_t i = 0; i < 20; i++) {
+ double dev = 0.5 * (inputBoundaries[i] + inputBoundaries[i + 1]) - 50;
+ inputData[i] =
+ exp(-dev * dev / 100) * (inputBoundaries[i + 1] - inputBoundaries[i]);
+ }
+ int indOfMax = VectorHelper::getBinIndex(inputBoundaries, 50.);
+ double fMax = inputData[indOfMax] /
+ (inputBoundaries[indOfMax + 1] - inputBoundaries[indOfMax]);
+ double iLeft = inputData[indOfMax - 1] /
+ (inputBoundaries[indOfMax] - inputBoundaries[indOfMax - 1]);
+ double iRight = inputData[indOfMax + 1] / (inputBoundaries[indOfMax + 2] -
+ inputBoundaries[indOfMax + 1]);
+
+ TS_ASSERT(iLeft < fMax);
+ TS_ASSERT(iRight < fMax);
+ VectorHelper::smoothInRange(inputData, output, 10, &inputBoundaries);
+ fMax = output[indOfMax] /
+ (inputBoundaries[indOfMax + 1] - inputBoundaries[indOfMax]);
+ iLeft = inputData[indOfMax - 1] /
+ (inputBoundaries[indOfMax] - inputBoundaries[indOfMax - 1]);
+ iRight = inputData[indOfMax + 1] /
+ (inputBoundaries[indOfMax + 2] - inputBoundaries[indOfMax + 1]);
+
+ TS_ASSERT(iLeft < fMax);
+ TS_ASSERT(iRight < fMax);
+
+ output.swap(inputData);
+ VectorHelper::smoothInRange(inputData, output, 10, &inputBoundaries);
+
+ fMax = output[indOfMax] /
+ (inputBoundaries[indOfMax + 1] - inputBoundaries[indOfMax]);
+ iLeft = inputData[indOfMax - 1] /
+ (inputBoundaries[indOfMax] - inputBoundaries[indOfMax - 1]);
+ iRight = inputData[indOfMax + 1] /
+ (inputBoundaries[indOfMax + 2] - inputBoundaries[indOfMax + 1]);
+
+ // TS_ASSERT(iLeft<fMax);
+ TS_ASSERT(iRight < fMax);
+
+ output.swap(inputData);
+ VectorHelper::smoothInRange(inputData, output, 10, &inputBoundaries);
+
+ fMax = output[indOfMax] /
+ (inputBoundaries[indOfMax + 1] - inputBoundaries[indOfMax]);
+ iLeft = inputData[indOfMax - 1] /
+ (inputBoundaries[indOfMax] - inputBoundaries[indOfMax - 1]);
+ iRight = inputData[indOfMax + 1] /
+ (inputBoundaries[indOfMax + 2] - inputBoundaries[indOfMax + 1]);
+
+ // TS_ASSERT(iLeft<fMax);
+ TS_ASSERT(iRight < fMax);
+
+ output.swap(inputData);
+ VectorHelper::smoothInRange(inputData, output, 10, &inputBoundaries);
+
+ fMax = output[indOfMax] /
+ (inputBoundaries[indOfMax + 1] - inputBoundaries[indOfMax]);
+ iLeft = inputData[indOfMax - 1] /
+ (inputBoundaries[indOfMax] - inputBoundaries[indOfMax - 1]);
+ iRight = inputData[indOfMax + 1] /
+ (inputBoundaries[indOfMax + 2] - inputBoundaries[indOfMax + 1]);
+
+ TS_ASSERT(inputData[indOfMax - 1] < output[indOfMax]);
+ TS_ASSERT(inputData[indOfMax + 1] < output[indOfMax]);
+ }
private:
/// Testing bins
diff --git a/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h b/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h
index c5d103d323ea76a43b145d302051846e0c8fcbcf..c2bbf05191f2bacbb0af8f49ad7d46f61e8d42bd 100644
--- a/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h
+++ b/Framework/TestHelpers/inc/MantidTestHelpers/WorkspaceCreationHelper.h
@@ -46,7 +46,7 @@ public:
return out;
}
};
-/** mock algorithn for doing logging/progress reporting*/
+/** mock algorithm for doing logging/progress reporting*/
class MockAlgorithm : public Mantid::API::Algorithm {
public:
MockAlgorithm(size_t nSteps = 100);
@@ -178,8 +178,8 @@ void addNoise(Mantid::API::MatrixWorkspace_sptr ws, double noise,
/**
* Create a test workspace with a fully defined instrument
* Each spectra will have a cylindrical detector defined 2*cylinder_radius away
- * from the centre of the
- * pervious.
+ * from the centre of the previous.
+ *
* Data filled with: Y: 2.0, E: sqrt(2.0), X: nbins of width 1 starting at 0
*/
Mantid::DataObjects::Workspace2D_sptr create2DWorkspaceWithFullInstrument(
@@ -262,7 +262,7 @@ CreateRandomEventWorkspace(size_t numbins, size_t numpixels,
Mantid::API::MatrixWorkspace_sptr
CreateGroupedWorkspace2D(size_t numHist, int numBins, double binDelta);
-// grouped workpsace with detectors arranges in rings in centre and into boxes
+// grouped workspace with detectors arranges in rings in center and into boxes
// outside
Mantid::API::MatrixWorkspace_sptr CreateGroupedWorkspace2DWithRingsAndBoxes(
size_t RootOfNumHist = 10, int numBins = 10, double binDelta = 1.0);
diff --git a/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp b/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp
index 88c67d797f42b4ec54f1238219a455cb436b9bf2..3d601ae49bc433f40bfb0c82aa8d02d8b936706c 100644
--- a/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp
+++ b/Framework/TestHelpers/src/WorkspaceCreationHelper.cpp
@@ -23,6 +23,7 @@
#include "MantidGeometry/Instrument/Detector.h"
#include "MantidGeometry/Instrument/ParameterMap.h"
#include "MantidGeometry/Instrument/ReferenceFrame.h"
+#include "MantidGeometry/Instrument/Component.h"
#include "MantidGeometry/Objects/ShapeFactory.h"
#include "MantidGeometry/Crystal/OrientedLattice.h"
#include "MantidKernel/MersenneTwister.h"
@@ -403,6 +404,10 @@ create2DWorkspaceWithFullInstrument(int nhist, int nbins, bool includeMonitors,
testInst->setPos(0.0, 0.0, 0.0);
testInst->add(sample);
testInst->markAsSamplePos(sample);
+ // chopper position
+ Component *chop_pos =
+ new Component("chopper-position", Kernel::V3D(-10, 0, 0), testInst.get());
+ testInst->add(chop_pos);
return space;
}
diff --git a/Testing/SystemTests/tests/analysis/ISISDirectReductionComponents.py b/Testing/SystemTests/tests/analysis/ISISDirectReductionComponents.py
index 214f517a7cc9a5b79b8f123ede6d7174597f4e39..437c4668bf1496371331dfb39d53947f1a15ef43 100644
--- a/Testing/SystemTests/tests/analysis/ISISDirectReductionComponents.py
+++ b/Testing/SystemTests/tests/analysis/ISISDirectReductionComponents.py
@@ -53,6 +53,7 @@ class ISIS_ReductionWebLike(stresstesting.MantidStressTest):
if 'outWS' in mtd:
return 'outWS'
saveFileName = self.rd.reducer.save_file_name
+#pylint: disable=unused-variable
outWS = Load(Filename=saveFileName+'.nxs')
outWS *= 0.997979227566217
fullRezPath =FileFinder.getFullPath(saveFileName+'.nxs')
@@ -63,7 +64,11 @@ class ISIS_ReductionWebLike(stresstesting.MantidStressTest):
def validate(self):
"""Returns the name of the workspace & file to compare"""
+ # tolerance defined outside of init
+#pylint: disable=W0201
self.tolerance = 1e-6
+ # tolerance_is_reller defined outside of init
+#pylint: disable=W0201
self.tolerance_is_reller=True
self.disableChecking.append('SpectraMap')
self.disableChecking.append('Instrument')
@@ -97,6 +102,8 @@ class ISIS_ReductionWrapperValidate(stresstesting.MantidStressTest):
# this is correct workflow for the ref file
#rd.reducer.prop_man.save_file_name = ref_file
# temporary workflow, until we fix workspace adjustment
+ # disable pylint -- access to protected member
+#pylint: disable=W0212
rd._tolerr =3.e-3
rd.reducer.prop_man.save_file_name = 'MARIReduction.nxs'
rd.validate_run_number=11001
@@ -239,6 +246,13 @@ class ISISLoadFilesMER(stresstesting.MantidStressTest):
det = mon_ws.getDetector(0)
self.assertTrue(det.isMonitor())
+ ei_ws = GetAllEi(mon_ws,69634,69638,IgnoreSecondMonitor=False)
+ self.assertTrue(isinstance(ei_ws,Workspace))
+
+ en_peaks = ei_ws.readX(0)
+ self.assertAlmostEquals(len(en_peaks),1)
+ self.assertAlmostEqual(en_peaks[0],108.94,2)
+
self.valid = True
@@ -283,12 +297,16 @@ class ISISLoadFilesLET(stresstesting.MantidStressTest):
#
self.assertEqual(mon_ws.getNumberHistograms(),27)
+ ei_ws = GetAllEi(mon_ws,40966,40967,IgnoreSecondMonitor=True)
+ self.assertTrue(isinstance(ei_ws,Workspace))
self.valid = True
+
def validate(self):
return self.valid
if __name__=="__main__":
- ISISLoadFilesMER.runTest()
+ tester = ISISLoadFilesLET()
+ tester.runTest()
diff --git a/Testing/SystemTests/tests/analysis/ValidateInstrumentDefinitionFiles.py b/Testing/SystemTests/tests/analysis/ValidateInstrumentDefinitionFiles.py
index f4df6a0b791a56ac268dde824c99b4884f2bf9e0..70f51e4427d150b0679b9fa6579287cf043b1e07 100644
--- a/Testing/SystemTests/tests/analysis/ValidateInstrumentDefinitionFiles.py
+++ b/Testing/SystemTests/tests/analysis/ValidateInstrumentDefinitionFiles.py
@@ -10,6 +10,8 @@ EXPECTED_EXT = '.expected'
class ValidateInstrumentDefinitionFiles(stresstesting.MantidStressTest):
xsdFile=''
+ # Explicitly specify single file to test. If None, test all.
+ theFileToTest=None #"MARI_Definition.xml"
def skipTests(self):
try:
@@ -66,7 +68,10 @@ class ValidateInstrumentDefinitionFiles(stresstesting.MantidStressTest):
direc = config['instrumentDefinition.directory']
self.xsdFile = os.path.join(direc,'Schema/IDF/1.0/','IDFSchema.xsd')
- files = self.__getDataFileList__()
+ if self.theFileToTest is None:
+ files = self.__getDataFileList__()
+ else:
+ files = [os.path.join(direc,self.theFileToTest)]
# run the tests
failed = []
@@ -94,5 +99,7 @@ class ValidateInstrumentDefinitionFiles(stresstesting.MantidStressTest):
if __name__ == '__main__':
valid = ValidateInstrumentDefinitionFiles()
+ # validate specific file
+ #valid.theFileToTest = "MARI_Definition.xml"
valid.runTest()
diff --git a/docs/source/algorithms/GetAllEi-v1.rst b/docs/source/algorithms/GetAllEi-v1.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c0e46da75eb0b368df50457ff64c7b2b4f1fb2e5
--- /dev/null
+++ b/docs/source/algorithms/GetAllEi-v1.rst
@@ -0,0 +1,95 @@
+.. algorithm::
+
+.. summary::
+
+.. alias::
+
+.. properties::
+
+Description
+-----------
+
+Algorithm finds the estimate for all incident energies allowed by chopper system of an inelastic instrument and returns a workspace,
+with the estimates for positions, heights and width of incident energies provided by the choppers and registered on monitors.
+These estimates can be used as guess value for :ref:`algm-GetEi` algorithm or as inputs for a peak fitting procedure.
+
+Algorithm performs number of steps to identify the values requested:
+
+#. It takes appropriate log names from instrument definition file (IDF), namely chopper-position component and calculates last chopper speed and delay as average
+ of the filtered log values. Guess chopper opening times are calculated from chopper speed and delay time. The "chopper-position" component with appropriate properties
+ has to be present in IDF for this algorithm to work. See ISIS MARI or MAPS instrument definition files for example of "chopper-position" component.
+
+#. Algorithm uses estimate for the minimal energy resolution of an instrument and searches for real peaks around guess values obtained
+ earlier within 4 sigma of this resolution interval.
+
+#. If peaks are found, the algorithm performs running averages over signal in the appropriate time interval until first derivative
+ of the signal has only one zero. The position of this zero is returned as the guess energy and the distance between closest to
+ the guess energy zeros of the second derivative are returned as the guess values for the peak width. The peak amplitude
+ is estimated from the total intensity of the signal within the search interval, assuming that the peak shape is Gaussian.
+
+#. Similar procedure is performed for second monitor. The peak is accepted only if the peak width lies between the minimal and maximal instrument resolution values
+ and the distance between peaks positions on two monitors (on energy scale) is smaller then two sigma.
+
+Algorithm returns matrix workspace containing single spectrum, with x-values representing peak positions, y-values: peak heights and the error: peak width. X-values are
+sorted according to energy in peaks (peaks with maximal energy are returned first).
+
+Used Subalgorithms
+------------------
+The algorithm uses :ref:`Unit Factory <Unit Factory>` and :ref:`algm-ConvertUnits` algorithm
+to convert units from TOF to energy.
+
+
+**Example: Find all incident energies for test workspace**
+
+.. testcode:: foundAllEi
+
+ # BUILD SAMPLE WORKSPACE
+ # Build sample workspace with chopper and in energy units to
+ # have defined peaks in defined energy positions
+ wsEn=CreateSampleWorkspace(Function='Multiple Peaks', NumBanks=1, BankPixelWidth=2, NumEvents=10000, XUnit='Energy', XMin=10, XMax=200, BinWidth=0.1)
+ # convert units to TOF to simulate real workspace obtained from experiment
+ ws = ConvertUnits(InputWorkspace=wsEn, Target='TOF')
+ # find chopper log values would be present in real workspace
+ l_chop = 7.5 # chopper position build into test workspace
+ l_mon1 = 15. # monitor 1 position (detector 1), build into test workspace
+ t_mon1 = 3100. # the time of flight defined by incident energy of the peak generated by CreateSampelpWorkspace algorithm.
+ t_chop = (l_chop/l_mon1)*t_mon1
+ # Add these log values to simulated workspace to represent real sample logs
+ AddTimeSeriesLog(ws, Name="fermi_delay", Time="2010-01-01T00:00:00", Value=t_chop ,DeleteExisting=True)
+ AddTimeSeriesLog(ws, Name="fermi_delay", Time="2010-01-01T00:30:00", Value=t_chop )
+ AddTimeSeriesLog(ws, Name="fermi_speed", Time="2010-01-01T00:00:00", Value=900 ,DeleteExisting=True)
+ AddTimeSeriesLog(ws, Name="fermi_speed", Time="2010-01-01T00:30:00", Value=900)
+ #-------------------------------------------------------------
+
+ # FIND GUESS PEAKS
+ allEiWs=GetAllEi(ws,Monitor1SpecID=1,Monitor2SpecID=2)
+ # Analyze results
+ allEi = allEiWs.readX(0);
+ peakHeight = allEiWs.readY(0);
+ peakWidth = allEiWs.readE(0);
+
+ # Check if peaks positions are indeed correct:
+ #-------------------------------------------------------------
+ resEi=[]
+ for ei_guess in allEi:
+ nop,t_peak,monIndex,tZero=GetEi(InputWorkspace=ws, Monitor1Spec=1, Monitor2Spec=2, EnergyEstimate=ei_guess)
+ resEi.append((nop,t_peak));
+ print "! Guess Ei ! peak TOF ! peak height ! peak width !"
+ for ind,val in enumerate(resEi):
+ print "! {0: >6.1f} ! {1: >6.2f} ! {2: >6.2f} ! {3: >6.2f} !".format(allEi[ind],val[1],peakHeight[ind],peakWidth[ind])
+ #
+ # NOTE: incident energy of GetEi is calculated from distance between monitor 1 and 2, and this distance is not correct in
+ # the test workspace. The tested point is that getEi can find energies from guess values and TOF for peaks is correct.
+
+Output:
+
+.. testoutput:: foundAllEi
+ :options: +NORMALIZE_WHITESPACE
+
+ ! Guess Ei ! peak TOF ! peak height ! peak width !
+ ! 67.0 ! 4188.03 ! 34.68 ! 2.35 !
+ ! 124.1 ! 3079.09 ! 14.01 ! 4.35 !
+
+.. categories::
+
+.. sourcelink::
diff --git a/instrument/LET_Definition.xml b/instrument/LET_Definition.xml
index a79bb314a74bb19549a4f00e59d65464e242a800..7080cc17821f2160e36408f35283559141c0a7f8 100644
--- a/instrument/LET_Definition.xml
+++ b/instrument/LET_Definition.xml
@@ -6,7 +6,7 @@
xsi:schemaLocation="http://www.mantidproject.org/IDF/1.0 http://schema.mantidproject.org/IDF/1.0/IDFSchema.xsd"
name="LET" valid-from ="1900-01-31 23:59:59"
valid-to ="2014-02-10 23:59:59"
- last-modified="2015-03-09 00:00:00">
+ last-modified="2015-03-09 00:00:00">
<defaults>
<length unit="meter"/>
@@ -672,6 +672,32 @@
<value units="microseconds" val="0"/>
</parameter>
</component-link>
+<!-- Chopper position -->
+ <component type="chopper-position">
+ <location z="-1.5"/>
+ <!--description is="The component provides sample-choper distance, used
+ to estimate chopper delay time and Tobyfit resolution calculations."/-->
+ <parameter name="initial_phase">
+ <value val="-0."/>
+ <description is="The initial rotation phase of the disk used to caluclate the time
+ for neutrons arriving at the chopper according to the formula time = delay + initial_phase/Speed"/>
+ </parameter>
+ <parameter name="ChopperDelayLog" type="string">
+ <value val="Chopper5_Disk2_phase"/>
+ </parameter>
+ <parameter name="ChopperSpeedLog" type="string">
+ <value val="Chopper5_Disk2_speed"/>
+ </parameter>
+ <parameter name="FilterBaseLog" type="string">
+ <value val="good_uah_log"/>
+ </parameter>
+ <!-- if the log above should be used as it is or
+ one should calculate its derivative -->
+ <parameter name="filter_with_derivative" type="bool">
+ <value val="True"/>
+ </parameter>
+ </component>
+ <type name="chopper-position" is="ChopperPos"></type>
<!-- Set the same across the rest of the instrument -->
<component-link name = "LET">
diff --git a/instrument/LET_Definition_dr2to12.xml b/instrument/LET_Definition_dr2to12.xml
index b3aaebc9f94ba4bf49a735c4af897b356ee0477a..063d32140eec5dc1f503d07252fb7045d2f94b1c 100644
--- a/instrument/LET_Definition_dr2to12.xml
+++ b/instrument/LET_Definition_dr2to12.xml
@@ -25,6 +25,7 @@
components defined in this file to face a position by default -->
<components-are-facing x="0.0" y="0.0" z="0.0" />
</defaults>
+
<!-- LIST OF PHYSICAL COMPONENTS (which the instrument consists of) -->
<!-- detector components -->
<properties>
@@ -89,6 +90,7 @@
</cylinder>
<algebra val="some-shape" />
</type>
+
<!-- Sample-position types -->
<type name="nickel-holder" is="SamplePos">
@@ -666,6 +668,35 @@
<id start="12470001" end="12471024" />
<id start="12480001" end="12481024" />
</idlist>
+<!-- Chopper position -->
+ <component type="chopper-position">
+ <location z="-1.5"/>
+ <!--description is="The component provides sample-choper distance, used
+ to estimate chopper delay time and Tobyfit resolution calculations."/-->
+ <parameter name="initial_phase">
+ <value val="-0."/>
+ <description is="The initial rotation phase of the disk used to caluclate the time
+ for neutrons arriving at the chopper according to the formula time = delay + initial_phase/Speed"/>
+ </parameter>
+ <parameter name="ChopperDelayLog" type="string">
+ <value val="Chopper5_Disk2_phase"/>
+ </parameter>
+ <parameter name="ChopperSpeedLog" type="string">
+ <value val="Chopper5_Disk2_speed"/>
+ </parameter>
+ <parameter name="FilterBaseLog" type="string">
+ <value val="good_uah_log"/>
+ </parameter>
+ <!-- if the log above should be used as it is or
+ one should calculate its derivative -->
+ <parameter name="filter_with_derivative" type="bool">
+ <value val="True"/>
+ </parameter>
+
+ </component>
+ <type name="chopper-position" is="ChopperPos"></type>
+
+
<!-- DETECTOR PARAMETERS -->
<component-link name="monitors">
<parameter name="DelayTime">
diff --git a/instrument/MAPS_Definition.xml b/instrument/MAPS_Definition.xml
index 5945f1ff1bdb314b8188dd48326e052f114087e9..ac7e3e16db7e15e5001ef00ef8a6f2cdbefc77c4 100644
--- a/instrument/MAPS_Definition.xml
+++ b/instrument/MAPS_Definition.xml
@@ -40,6 +40,36 @@
<location />
</component>
<type name="sample-position" is="SamplePos"></type>
+
+<!-- Chopper position -->
+ <component type="chopper-position">
+ <location z="-1.895"/>
+ <!--description is="The component provides sample-choper distance, used
+ to estimate chopper delay time and Tobyfit resolution calculations."/-->
+ <parameter name="initial_phase">
+ <value val="130000."/>
+ <description is="The initial rotation phase of the disk used to caluclate the time
+ for neutrons arriving at the chopper according to the formula time = delay + initial_phase/Speed"/>
+ </parameter>
+ <parameter name="ChopperDelayLog" type="string">
+ <value val="Chopper_delay"/>
+ </parameter>
+ <parameter name="ChopperSpeedLog" type="string">
+ <value val="Chopper_Speed"/>
+ </parameter>
+ <parameter name="FilterBaseLog" type="string">
+ <value val="good_uah_log"/>
+ </parameter>
+ <!-- if the log above should be used as it is or
+ one should calculate its derivative -->
+ <parameter name="filter_with_derivative" type="bool">
+ <value val="True"/>
+ </parameter>
+
+ </component>
+ <type name="chopper-position" is="ChopperPos"></type>
+
+
<component type="monitors" idlist="monitors">
diff --git a/instrument/MARI_Definition.xml b/instrument/MARI_Definition.xml
index 41a9900ffe0290b5c789aa22b3bd8cc6bd0e2d95..5ef895cb7b30160d062f8001a44364d251267eb1 100644
--- a/instrument/MARI_Definition.xml
+++ b/instrument/MARI_Definition.xml
@@ -42,6 +42,38 @@
<!-- CHOPPERS -->
+ <!-- Chopper position -->
+ <component type="chopper-position">
+ <location z="-1.689"/>
+ <!--description is="The component provides sample-choper distance, used
+ to estimate chopper delay time and Tobyfit resolution calculations."/ -->
+ <!-- if chopper "a": t_offset = -3281.8/SPEED + 16.2;
+ if chopper "g": t_offset = -2822.5/SPEED + 9.7
+ if chopper "s": t_offset = -3326.8/SPEED + 7.7
+ -->
+ <parameter name="initial_phase">
+ <value val="-3000."/>
+ <description is="The initial rotation phase of the disk used to caluclate the time
+ for neutrons arriving at the chopper according to the formula time = delay + initial_phase/Speed"/>
+ </parameter>
+ <parameter name="ChopperDelayLog" type="string">
+ <value val="fermi_delay"/>
+ </parameter>
+ <parameter name="ChopperSpeedLog" type="string">
+ <value val="fermi_speed"/>
+ </parameter>
+ <parameter name="FilterBaseLog" type="string">
+ <value val="good_uah_log"/>
+ </parameter>
+ <!-- if the log above should be used as it is or
+ one should calculate its derivative -->
+ <parameter name="filter_with_derivative" type="bool">
+ <value val="True"/>
+ </parameter>
+
+ </component>
+ <type name="chopper-position" is="ChopperPos"></type>
+
<component type="fermi-chopper">
<location y="-0.1" z="-1.689" />
@@ -50,9 +82,11 @@
</parameter>
<parameter name="Delay (us)">
<logfile id="fermi_delay" extract-single-value-as="last_value" />
- <description is="Delay is the time chopper waits for...? (test description)"/>
+ <description is="Delay is the time between the moment when beam hits moderator and the moment,
+ when neutron passes through chopper for the first time."/>
</parameter>
- </component>
+
+ </component>
<type name="fermi-chopper" is="chopper">
<cylinder id="body">
<centre-of-bottom-base x="0.0" y="0.0" z="0.0" />
diff --git a/instrument/MERLIN_Definition_after2013_4.xml b/instrument/MERLIN_Definition_after2013_4.xml
index fd69107e5c7fdd2600996cd426f5dfb5fb84af57..4bb893e054a38bf4804f602fffac62d6dc918b69 100644
--- a/instrument/MERLIN_Definition_after2013_4.xml
+++ b/instrument/MERLIN_Definition_after2013_4.xml
@@ -62,8 +62,31 @@
<!-- Chopper position -->
<component type="chopper-position">
<location z="-1.8"/>
+ <!--description is="The component provides sample-choper distance, used
+ to estimate chopper delay time and Tobyfit resolution calculations."/-->
+ <parameter name="initial_phase">
+ <value val="-490000."/>
+ <description is="The initial rotation phase of the disk used to caluclate the time
+ for neutrons arriving at the chopper according to the formula time = delay + initial_phase/Speed"/>
+ </parameter>
+ <parameter name="ChopperDelayLog" type="string">
+ <value val="Chopper_delay"/>
+ </parameter>
+ <parameter name="ChopperSpeedLog" type="string">
+ <value val="Chopper_Speed"/>
+ </parameter>
+ <parameter name="FilterBaseLog" type="string">
+ <value val="good_uah_log"/>
+ </parameter>
+ <!-- if the log above should be used as it is or
+ one should calculate its derivative -->
+ <parameter name="filter_with_derivative" type="bool">
+ <value val="True"/>
+ </parameter>
+
</component>
- <type name="chopper-position" is="ChopperPos"></type>
+ <type name="chopper-position" is="ChopperPos"></type>
+
<!-- DETECTORS -->
diff --git a/instrument/Schema/IDF/1.0/IDFSchema.xsd b/instrument/Schema/IDF/1.0/IDFSchema.xsd
index 4b2ea38bf59fbf02f4d836c6ae2b28c8e0c6e0f5..d2a3c79d955145c7a7aeef77f51e1ef3634fff75 100644
--- a/instrument/Schema/IDF/1.0/IDFSchema.xsd
+++ b/instrument/Schema/IDF/1.0/IDFSchema.xsd
@@ -289,9 +289,11 @@
<xs:simpleContent>
<xs:extension base="xs:string">
<xs:anyAttribute/>
+
</xs:extension>
</xs:simpleContent>
</xs:complexType>
+
</xs:element>
<xs:element name="properties">
<xs:complexType>