diff --git a/Code/Mantid/Algorithms/inc/MantidAlgorithms/GetEi2.h b/Code/Mantid/Algorithms/inc/MantidAlgorithms/GetEi2.h
index 0f438c9a5e4556ca4c8d79131ccea308278cbc88..552618e8c48eb4d88209ed28776c92a6184101c2 100644
--- a/Code/Mantid/Algorithms/inc/MantidAlgorithms/GetEi2.h
+++ b/Code/Mantid/Algorithms/inc/MantidAlgorithms/GetEi2.h
@@ -5,11 +5,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/Algorithm.h"
-#include "MantidDataObjects/Workspace2D.h"
-#include "MantidGeometry/Instrument/ParametrizedComponent.h"
-#include "MantidGeometry/Instrument/Component.h"
-#include <vector>
-#include <string>
+#include "MantidAPI/MatrixWorkspace.h"
namespace Mantid
{
@@ -53,9 +49,8 @@ namespace Algorithms
class DLLExport GetEi2 : public API::Algorithm
{
public:
- /// empty contructor calls the base class constructor
- GetEi2() : Algorithm(), m_mon_indices() {}
-
+ /// Default constructor
+ GetEi2();
/// Initialize the algorithm
void init();
/// Execute the algorithm
@@ -69,30 +64,48 @@ namespace Algorithms
virtual const std::string category() const{return "CorrectionFunctions";}
private:
- void advanceProgress(double toAdd);
- void getGeometry(DataObjects::Workspace2D_const_sptr WS, double &monitor0Dist, double &monitor1Dist) const;
- double timeToFly(double s, double E_KE) const;
- double getPeakCentre(API::MatrixWorkspace_const_sptr WS, const int monitIn, const double peakTime);
- void extractSpec(int specInd, double start = -1.0, double end = -1.0);
- double getPeakFirstMoments(API::MatrixWorkspace_sptr WS, const double tMin, const double tMax);
- void regroup(double xmin, double delta, double xmax, const MantidVec &x, const MantidVec &y, const MantidVec &e, MantidVec& xnew, MantidVec& ynew, MantidVec& enew);
- void getPeakMean(const MantidVec& Xs, const MantidVec& Ys, const MantidVec& Es, const double prominence, double &area, double &c, double &c_fwhm, double &w, double &xbar);
- void integrate(double &bkgd_m, double &bkgd_err_m, const MantidVec &x, const MantidVec &y, const MantidVec &e, const double start, const double end);
- API::MatrixWorkspace_sptr smooth(API::MatrixWorkspace_sptr WS);
- API::MatrixWorkspace_sptr rebin(API::MatrixWorkspace_sptr WS, const double first, const double width, const double end);
- double neutron_E_At(double speed) const;
- /// An estimate of the percentage of the algorithm runtimes that has been completed
- double m_fracCompl;
- /// name of the tempory workspace that we create and will contain the monitor histogram that we're examining
- API::MatrixWorkspace_sptr m_tempWS;
- /// Workspace indices for the monitors
- std::vector<int> m_mon_indices;
-
- // for estimating algorithm progress
- static const double CROP; ///< fraction of algorithm time taken up with running CropWorkspace
- static const double GET_COUNT_RATE; ///< fraction of algorithm taken by a single call to ConvertToDistribution
- static const double FIT_PEAK; ///< fraction required to find a peak
+ /// Calculate Ei from the initial guess given
+ double calculateEi(const double initial_guess);
+ /// Get the distance from the source of the detector at the workspace index given
+ double getDistanceFromSource(const int ws_index) const;
+ /// Calculate the peak position within the given window
+ double calculatePeakPosition(const int ws_index, const double t_min, const double t_max);
+ /// Extract the required spectrum from the input workspace
+ API::MatrixWorkspace_sptr extractSpectrum(const int ws_index, const double start, const double end);
+ /// Calculate peak width
+ double calculatePeakWidthAtHalfHeight(API::MatrixWorkspace_sptr data_ws, const double prominence,
+ std::vector<double> & peak_x, std::vector<double> & peak_y, std::vector<double> & peak_e) const;
+ /// Calculate the value of the first moment of the given spectrum
+ double calculateFirstMoment(API::MatrixWorkspace_sptr monitor_ws, const double prominence);
+ /// Rebin the given workspace using the given parameters
+ API::MatrixWorkspace_sptr rebin(API::MatrixWorkspace_sptr input_ws, const double first, const double width, const double end);
+ /// Integrate the point data
+ void integrate(double &integral_value, double &integral_err, const MantidVec &x, const MantidVec &y, const MantidVec &e, const double xmin, const double xmax) const;
+ /// Offset the bins on the input workspace by the calculated time of the first monitor peak
+ void applyBinOffset();
+ /// Move the source of neutrons to the position defined by the first input monitor
+ void moveNeutronSource();
+ /// Store the incident energy within the sample object
+ void storeEi(const double ei) const;
+ /// The input workspace
+ API::MatrixWorkspace_sptr m_input_ws;
+ /// The calculated position of the first peak
+ std::pair<int, double> m_peak1_pos;
+ /// True if the Ei should be fixed at the guess energy
+ bool m_fixedei;
+ /// Conversion factor between time and energy
+ double m_t_to_mev;
+ /// The percentage deviation from the estimated peak time that defines the peak region
+ const double m_tof_window;
+ /// Number of std deviations to consider a peak
+ const double m_peak_signif;
+ /// Number of std deviations to consider a peak for the derivative
+ const double m_peak_deriv;
+ /// The fraction of the peak width for the new bins
+ const double m_binwidth_frac;
+ /// The fraction of the peakwidth to use in calculating background points
+ const double m_bkgd_frac;
};
} // namespace Algorithms
diff --git a/Code/Mantid/Algorithms/src/GetEi2.cpp b/Code/Mantid/Algorithms/src/GetEi2.cpp
index 650fbdb763d038bc8a369e94816ce0a250ce37f0..0034ded846c92d2a6efd7f55c04d76564506484f 100644
--- a/Code/Mantid/Algorithms/src/GetEi2.cpp
+++ b/Code/Mantid/Algorithms/src/GetEi2.cpp
@@ -1,688 +1,689 @@
#include "MantidAlgorithms/GetEi2.h"
-#include "MantidKernel/ArrayProperty.h"
-#include "MantidKernel/FileProperty.h"
+
#include "MantidKernel/PhysicalConstants.h"
+#include "MantidKernel/VectorHelper.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidGeometry/Instrument/DetectorGroup.h"
-#include "MantidKernel/Exception.h"
-#include "MantidKernel/VectorHelper.h"
+#include "MantidGeometry/IObjComponent.h"
+#include "MantidDataObjects/Workspace2D.h"
#include <boost/lexical_cast.hpp>
-
#include <cmath>
#include <algorithm>
-#include <iostream>
#include <sstream>
-#include <iomanip>
namespace Mantid
{
namespace Algorithms
{
-
// Register the algorithm into the algorithm factory
DECLARE_ALGORITHM(GetEi2)
-
- using namespace Kernel;
- using namespace API;
- using namespace Geometry;
- using namespace DataObjects;
-
- // @todo: Move to statics on the class
- static const double MON_TOF_WIN=0.1;
- static const double SFAC_PEAK = 2.0, SFAC_DERIV=1.0;
- static const int NPOINTS=8;
- static const double BKGD_FAC=0.5;
- static const double NO_ESTIMATE=-1e200;
-
- // progress estimates
- const double GetEi2::CROP = 0.15;
- const double GetEi2::GET_COUNT_RATE = 0.15;
- const double GetEi2::FIT_PEAK = 0.2;
-
-
- void GetEi2::init()
-
- {// Declare required input parameters for algorithm and do some validation here
- CompositeValidator<Workspace2D> *val = new CompositeValidator<Workspace2D>;
- val->add(new WorkspaceUnitValidator<Workspace2D>("TOF"));
- val->add(new HistogramValidator<Workspace2D>);
- declareProperty(new WorkspaceProperty<Workspace2D>(
- "InputWorkspace","",Direction::Input,val),
- "The X units of this workspace must be time of flight with times in\n"
- "micro-seconds");
- BoundedValidator<int> *mustBePositive = new BoundedValidator<int>();
- mustBePositive->setLower(0);
- declareProperty("Monitor1Spec", -1, mustBePositive,
- "The spectrum number of the output of the first monitor, e.g. MAPS\n"
- "41474, MARI 2, MERLIN 69634");
- declareProperty("Monitor2Spec", -1, mustBePositive->clone(),
- "The spectrum number of the output of the second monitor e.g. MAPS\n"
- "41475, MARI 3, MERLIN 69638");
- BoundedValidator<double> *positiveDouble = new BoundedValidator<double>();
- positiveDouble->setLower(0.0);
- declareProperty("EnergyEstimate", EMPTY_DBL(), positiveDouble,
- "An approximate value for the typical incident energy, energy of\n"
- "neutrons leaving the source (meV)");
- declareProperty("IncidentEnergy", -1.0, Direction::Output);
- declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
-
- m_fracCompl = 0.0;
- }
-
- /** Executes the algorithm
- * @throw out_of_range if the peak runs off the edge of the histogram
- * @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
- * @throw IndexError if there is a problem converting spectra indexes to spectra numbers, which would imply there is a problem with the workspace
- * @throw invalid_argument if a good peak fit wasn't made or the input workspace does not have common binning
- * @throw runtime_error if there is a problem with the SpectraDetectorMap or a sub-algorithm falls over
- */
- void GetEi2::exec()
+ }
+}
+
+using namespace Mantid::Kernel;
+using namespace Mantid::API;
+using namespace Mantid::Geometry;
+using namespace Mantid::Algorithms;
+
+/**
+* Default contructor
+*/
+GetEi2::GetEi2() : Algorithm(), m_input_ws(), m_peak1_pos(0, 0.0), m_fixedei(false), m_tof_window(0.1), m_peak_signif(2.0), m_peak_deriv(1.0),
+ m_binwidth_frac(1.0/12.0), m_bkgd_frac(0.5)
+{
+ // Conversion factor common for converting between micro seconds and energy in meV
+ m_t_to_mev = 5e11 * PhysicalConstants::NeutronMass / PhysicalConstants::meV;
+}
+
+void GetEi2::init()
+
+{// Declare required input parameters for algorithm and do some validation here
+ CompositeValidator<MatrixWorkspace> *validator = new CompositeValidator<MatrixWorkspace>;
+ validator->add(new WorkspaceUnitValidator<MatrixWorkspace>("TOF"));
+ validator->add(new HistogramValidator<MatrixWorkspace>);
+
+ declareProperty(new WorkspaceProperty<MatrixWorkspace>(
+ "InputWorkspace","",Direction::InOut, validator),
+ "The X units of this workspace must be time of flight with times in\n"
+ "microseconds");
+ BoundedValidator<int> *mustBePositive = new BoundedValidator<int>();
+ mustBePositive->setLower(0);
+ declareProperty("Monitor1Spec", -1, mustBePositive,
+ "The spectrum number to use as the first monitor\n");
+ declareProperty("Monitor2Spec", -1, mustBePositive->clone(),
+ "The spectrum number to use as the second monitor\n");
+ BoundedValidator<double> *positiveDouble = new BoundedValidator<double>();
+ positiveDouble->setLower(0.0);
+ declareProperty("EnergyEstimate", -1.0 , positiveDouble,
+ "An approximate value for the typical incident energy, energy of\n"
+ "neutrons leaving the source (meV)");
+ declareProperty("FixEi", false, "If true, the incident energy will be set to the value of the \n"
+ "EnergyEstimate property.");
+ declareProperty("AdjustBins", false, "If true, the bins from the input workspace \n"
+ "will be offset by the first monitor peak and the source moved to the position of \n"
+ "the monitor defined by Monitor1Spec.");
+ declareProperty("IncidentEnergy", -1.0, Direction::Output);
+ declareProperty("FirstMonitorPeak", -1.0, Direction::Output);
+}
+
+/** Executes the algorithm
+* @throw out_of_range if the peak runs off the edge of the histogram
+* @throw NotFoundError if one of the requested spectrum numbers was not found in the workspace
+* @throw IndexError if there is a problem converting spectra indexes to spectra numbers, which would imply there is a problem with the workspace
+* @throw invalid_argument if a good peak fit wasn't made or the input workspace does not have common binning
+* @throw runtime_error if there is a problem with the SpectraDetectorMap or a sub-algorithm falls over
+*/
+void GetEi2::exec()
+{
+ m_input_ws = getProperty("InputWorkspace");
+ m_fixedei = getProperty("FixEi");
+ double initial_guess = getProperty("EnergyEstimate");
+ double incident_energy = calculateEi(initial_guess);
+ if( !m_fixedei )
+ {
+ g_log.notice() << "Incident energy = " << incident_energy << " meV from initial guess = " << initial_guess << " meV\n";
+ }
+
+ storeEi(incident_energy);
+
+ const bool adjust_bins = getProperty("AdjustBins");
+ if( adjust_bins )
+ {
+ try
{
- Workspace2D_const_sptr inWS = getProperty("InputWorkspace");
- const int monitor1_spec = getProperty("Monitor1Spec");
- const int monitor2_spec = getProperty("Monitor2Spec");
-
- //Covert spectrum numbers to workspace indices
- std::vector<int> spec_nums(2, monitor1_spec);
- spec_nums[1] = monitor2_spec;
- m_mon_indices.clear();
- // get the index number of the histogram for the first monitor
- WorkspaceHelpers::getIndicesFromSpectra(inWS, spec_nums, m_mon_indices);
-
- if( m_mon_indices.size() != 2 )
- {
- g_log.error() << "Error retrieving monitor spectra from input workspace. Check input properties.\n";
- throw std::runtime_error("Error retrieving monitor spectra spectra from input workspace.");
- }
-
- double dist2moni0(-1.0), dist2moni1(-1.0);
- getGeometry(inWS, dist2moni0, dist2moni1);
- g_log.debug() << "Distance between monitors = " << dist2moni0 - dist2moni1 << " m\n";
-
- // the E_i estimate is used to find (identify) the monitor peaks, checking prior to fitting will throw an exception if this estimate is too big or small
- const double E_est = getProperty("EnergyEstimate");
- double peakLoc0(0.0);
- if( E_est != EMPTY_DBL() )
- {
- peakLoc0 = 1e6*timeToFly(dist2moni0, E_est);
- }
-
- g_log.information() << "Looking for a peak in the first monitor spectrum, workspace index " << m_mon_indices[0] << std::endl;
- double t_monitor0 = getPeakCentre(inWS, 1, peakLoc0);
- g_log.notice() << "The first peak has been found at TOF = " << t_monitor0 << " microseconds\n";
- setProperty("FirstMonitorPeak", t_monitor0);
-
- g_log.information() << "Looking for a peak in the second monitor spectrum, workspace index " << m_mon_indices[1] << std::endl;
- const double peakLoc1 = 1e6*timeToFly(dist2moni1, E_est);
- double t_monitor1 = getPeakCentre(inWS, 2, peakLoc1);
- g_log.information() << "The second peak has been found at TOF = " << t_monitor1 << " microseconds\n";
-
- // assumes that the source and the both mintors lie on one straight line, the 1e-6 converts microseconds to seconds as the mean speed needs to be in m/s
- double meanSpeed = 1e6 * (dist2moni1 - dist2moni0)/ (t_monitor1 - t_monitor0);
-
- // uses 0.5mv^2 to get the kinetic energy in joules which we then convert to meV
- double E_i = neutron_E_At(meanSpeed)/PhysicalConstants::meV;
- g_log.notice() << "The incident energy has been calculated to be " << E_i << " meV";
- if( E_est != EMPTY_DBL() )
- {
- g_log.information () << " (your estimate was " << E_est << " meV)\n";
- }
- else
- {
- g_log.information () << " (No estimate was supplied).\n";
- }
- setProperty("IncidentEnergy", E_i);
- }
-
- /** Gets the distances between the source and detectors whose IDs you pass to it
- * @param WS the input workspace
- * @param monitor0Dist the calculated distance to the detector whose ID was passed to this function first
- * @param monitor1Dist calculated distance to the detector whose ID was passed to this function second
- * @throw NotFoundError if no detector is found for the detector ID given
- * @throw runtime_error if there is a problem with the SpectraDetectorMap
- */
- void GetEi2::getGeometry(DataObjects::Workspace2D_const_sptr WS, double &monitor0Dist, double &monitor1Dist) const
- {
- const IObjComponent_sptr source = WS->getInstrument()->getSource();
-
- // retrieve a pointer to the first monitor and get its distance
- IDetector_sptr det = WS->getDetector(m_mon_indices.front());
- if( !det )
- {
- g_log.error() << "A detector for monitor at workspace index " << m_mon_indices.front() << " cannot be found. ";
- throw std::runtime_error("No detector found for the first monitor.");
- }
- if( boost::dynamic_pointer_cast<DetectorGroup>(det) )
- {
- g_log.error() << "The detector for spectrum number " << m_mon_indices.front() << " is a group, grouped monitors are not supported by this algorithm\n";
- g_log.error() << "Error retrieving data for the first monitor" << std::endl;
- throw std::runtime_error("Detector for first monitor is a DetectorGroup.");
- }
- monitor0Dist = det->getDistance(*source);
- g_log.information() << "L1 distance " << monitor0Dist << " m.\n";
-
- // second monitor
- det = WS->getDetector(m_mon_indices[1]);
- if( !det )
- {
- g_log.error() << "A detector for monitor at workspace index " << m_mon_indices[1] << " cannot be found. ";
- throw std::runtime_error("No detector found for the second monitor.");
- }
- if( boost::dynamic_pointer_cast<DetectorGroup>(det) )
- {
- g_log.error() << "The detector for spectrum number " << m_mon_indices.front() << " is a group, grouped monitors are not supported by this algorithm\n";
- g_log.error() << "Error retrieving data for the second monitor" << std::endl;
- throw std::runtime_error("Detector for second monitor is a DetectorGroup.");
- }
- monitor1Dist = det->getDistance(*source);
- g_log.information() << "L2 distance " << monitor1Dist << " m.\n";
+ applyBinOffset();
+ moveNeutronSource();
}
-
- /** Uses E_KE = mv^2/2 and s = vt to calculate the time required for a neutron
- * to travel a distance, s
- * @param s ditance travelled in meters
- * @param E_KE kinetic energy in meV
- * @return the time to taken to travel that uninterrupted distance in seconds
- */
- double GetEi2::timeToFly(double s, double E_KE) const
- {
- // E_KE = mv^2/2, s = vt
- // t = s/v, v = sqrt(2*E_KE/m)
- // t = s/sqrt(2*E_KE/m)
-
- // convert E_KE to joules kg m^2 s^-2
- E_KE *= PhysicalConstants::meV;
-
- return s/sqrt(2*E_KE/PhysicalConstants::NeutronMass);
+ catch(std::exception& e)
+ {
+ g_log.warning() << "Bin offsetting failed. " << e.what() << "\n";;
}
-
- /** Looks for and examines a peak close to that specified by the input parameters and
- * examines it to find a representative time for when the neutrons hit the detector
- * @param WS the workspace containing the monitor spectrum
- * @param peak First or second monitor peak.
- * @param peakTime the estimated TOF of the monitor peak in the time units of the workspace
- * @return a time of flight value in the peak in microseconds
- * @throw invalid_argument if a good peak fit wasn't made or the input workspace does not have common binning
- * @throw out_of_range if the peak runs off the edge of the histogram
- * @throw runtime_error a sub-algorithm just falls over
- */
- double GetEi2::getPeakCentre(API::MatrixWorkspace_const_sptr WS, const int peak, const double peakTime)
+ }
+
+ setProperty("InputWorkspace", m_input_ws);
+ setProperty("IncidentEnergy", incident_energy);
+ setProperty("FirstMonitorPeak", m_peak1_pos.second);
+}
+
+/**
+ * Calculate the incident energy of the neutrons for the input workspace on this algorithm
+ * @param initial_guess A guess for value of the incident energy
+ */
+double GetEi2::calculateEi(const double initial_guess)
+{
+ const int monitor1_spec = getProperty("Monitor1Spec");
+ const int monitor2_spec = getProperty("Monitor2Spec");
+
+ //Covert spectrum numbers to workspace indices
+ std::vector<int> spec_nums(2, monitor1_spec);
+ spec_nums[1] = monitor2_spec;
+ std::vector<int> mon_indices;
+ mon_indices.reserve(2);
+ // get the index number of the histogram for the first monitor
+ WorkspaceHelpers::getIndicesFromSpectra(m_input_ws, spec_nums, mon_indices);
+
+ if( mon_indices.size() != 2 )
+ {
+ g_log.error() << "Error retrieving monitor spectra from input workspace. Check input properties.\n";
+ throw std::runtime_error("Error retrieving monitor spectra spectra from input workspace.");
+ }
+
+ // Calculate actual peak postion for each monitor peak
+ double peak_times[2] = {0.0, 0.0};
+ double det_distances[2] = {0.0, 0.0};
+ const bool adjust_bins = getProperty("AdjustBins");
+ for( unsigned int i = 0; i < 2; ++i )
+ {
+ int ws_index = mon_indices[i];
+ det_distances[i] = getDistanceFromSource(ws_index);
+ const double peak_guess = det_distances[i]*std::sqrt(m_t_to_mev/initial_guess);
+ if( m_fixedei && i == 0 )
{
- double tMin(0.0), tMax(0.0);
-
- if ( peakTime > 0.0 )
- {
- tMin = (1 - MON_TOF_WIN)*peakTime;
- tMax = (1 + MON_TOF_WIN)*peakTime;
- switch (peak)
- {
- case 1:
- g_log.information() << "Based on the user selected energy the first peak will be searched for at TOF " << peakTime << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*MON_TOF_WIN) << "%\n";
- break;
- case 2:
- g_log.information() << "Based on the user selected energy the second peak will be searched for at TOF " << peakTime << " micro seconds +/-" << boost::lexical_cast<std::string>(100.0*MON_TOF_WIN) << "%\n";
- break;
- default:
- throw std::runtime_error("Invalid monitor selected");
- }
- }
- else
- {
- tMin = 400;
- tMax = 12000;
- g_log.information() << "No energy estimate given, using default window t0 = " << tMin << " microseconds, t1 = " << tMax << " microseconds\n";
- }
- // runs CropWorkspace as a sub-algorithm to and puts the result in a new temporary workspace that will be deleted when this algorithm has finished
- extractSpec(m_mon_indices[peak-1], tMin, tMax);
- // converting the workspace to count rate is required by the fitting algorithm if the bin widths are not all the same, if the workspace is already a distribution this does nothing
- WorkspaceHelpers::makeDistribution(m_tempWS);
-
- return getPeakFirstMoments(m_tempWS, tMin, tMax);
+ m_peak1_pos = std::make_pair(ws_index, peak_guess);
+ g_log.information() << "First monitor peak = " << peak_guess << " microseconds from fixed Ei = " << initial_guess << " meV\n";
+ break;
}
- /** Calls CropWorkspace as a sub-algorithm and passes to it the InputWorkspace property
- * @param specInd the index number of the histogram to extract
- * @param start the number of the first bin to include (starts counting bins at 0)
- * @param end the number of the last bin to include (starts counting bins at 0)
- * @throw out_of_range if start, end or specInd are set outside of the vaild range for the workspace
- * @throw runtime_error if the algorithm just falls over
- * @throw invalid_argument if the input workspace does not have common binning
- */
- void GetEi2::extractSpec(int specInd, double start, double end)
- {
- IAlgorithm_sptr childAlg =
- createSubAlgorithm("CropWorkspace", 100*m_fracCompl, 100*(m_fracCompl+CROP) );
- m_fracCompl += CROP;
-
- childAlg->setPropertyValue( "InputWorkspace",
- getPropertyValue("InputWorkspace") );
- childAlg->setProperty( "StartWorkspaceIndex", specInd);
- childAlg->setProperty( "EndWorkspaceIndex", specInd);
+ const double t_min = (1.0 - m_tof_window)*peak_guess;
+ const double t_max = (1.0 + m_tof_window)*peak_guess;
+ g_log.information() << "Time-of-flight window for peak " << (i+1) << ": tmin = " << t_min << " microseconds, tmax = " << t_max << "microseconds\n";
+ peak_times[i] = calculatePeakPosition(ws_index, t_min, t_max);
+ g_log.information() << "Peak for monitor " << (i+1) << " = " << peak_times[i] << " microseconds\n";
+ if( i == 0 )
+ {
+ //Store for later adjustment of bins
+ m_peak1_pos = std::make_pair(ws_index, peak_times[i]);
+ }
+ }
+
+ if( m_fixedei )
+ {
+ return initial_guess;
+ }
+ else
+ {
+ double mean_speed = (det_distances[1] - det_distances[0])/(peak_times[1] - peak_times[0]);
+ return mean_speed*mean_speed*m_t_to_mev;
+ }
+}
+
+/** Gets the distance between the source and detectors whose workspace index is passed
+* @param ws_index The workspace index of the detector
+* @throw runtime_error if there is a problem
+*/
+double GetEi2::getDistanceFromSource(int ws_index) const
+{
+ const IObjComponent_sptr source = m_input_ws->getInstrument()->getSource();
+ // Retrieve a pointer detector
+ IDetector_sptr det = m_input_ws->getDetector(ws_index);
+ if( !det )
+ {
+ g_log.error() << "A detector for monitor at workspace index " << ws_index << " cannot be found. ";
+ throw std::runtime_error("No detector found.");
+ }
+ if( boost::dynamic_pointer_cast<DetectorGroup>(det) )
+ {
+ g_log.error() << "The detector for workspace " << ws_index << " is a group, grouped monitors are not supported by this algorithm\n";
+ throw std::runtime_error("Detector for monitor is a DetectorGroup.");
+ }
+ return det->getDistance(*source);
+}
+
+/**
+ * Calculate the peak position of a given TOF range within the chosen spectrum
+ */
+double GetEi2::calculatePeakPosition(int ws_index, double t_min, double t_max)
+{
+ //Crop out the current monitor workspace to the min/max times defined
+ MatrixWorkspace_sptr monitor_ws = extractSpectrum(ws_index, t_min, t_max);
+ // Workspace needs to be a count rate for the fitting algorithm
+ WorkspaceHelpers::makeDistribution(monitor_ws);
- if( start > 0 )
- {
- childAlg->setProperty( "XMin", start);
- childAlg->setProperty( "XMax", end);
- }
+ const double prominence(4.0);
+ std::vector<double> dummyX, dummyY, dummyE;
+ double peak_width = calculatePeakWidthAtHalfHeight(monitor_ws, prominence, dummyX, dummyY, dummyE);
+ monitor_ws = rebin(monitor_ws, t_min, peak_width*m_binwidth_frac, t_max);
- try
- {
- childAlg->execute();
- }
- catch (std::exception&)
- {
- g_log.error("Exception thrown while running CropWorkspace as a sub-algorithm");
- throw;
- }
+ double t_mean(0.0);
+ try
+ {
+ t_mean = calculateFirstMoment(monitor_ws, prominence);
+ }
+ catch(std::runtime_error&)
+ {
+ //Retry with smaller prominence factor
+ t_mean = calculateFirstMoment(monitor_ws, 2.0);
+ }
- if ( ! childAlg->isExecuted() )
- {
- g_log.error("The CropWorkspace algorithm failed unexpectedly, aborting.");
- throw std::runtime_error(name() + " failed trying to run CropWorkspace");
- }
- m_tempWS = childAlg->getProperty("OutputWorkspace");
+ return t_mean;
+}
+
+/** Calls CropWorkspace as a sub-algorithm and passes to it the InputWorkspace property
+* @param ws_index the index number of the histogram to extract
+* @param start the number of the first bin to include (starts counting bins at 0)
+* @param end the number of the last bin to include (starts counting bins at 0)
+* @throw out_of_range if start, end or specInd are set outside of the vaild range for the workspace
+* @throw runtime_error if the algorithm just falls over
+* @throw invalid_argument if the input workspace does not have common binning
+*/
+MatrixWorkspace_sptr GetEi2::extractSpectrum(int ws_index, const double start, const double end)
+{
+ IAlgorithm_sptr childAlg = createSubAlgorithm("CropWorkspace");
+ DataObjects::Workspace2D_sptr input_2D = boost::dynamic_pointer_cast<DataObjects::Workspace2D>(m_input_ws);
+ childAlg->setProperty("InputWorkspace", input_2D);
+ childAlg->setProperty("StartWorkspaceIndex", ws_index);
+ childAlg->setProperty("EndWorkspaceIndex", ws_index);
+ childAlg->setProperty("XMin", start);
+ childAlg->setProperty("XMax", end);
+ try
+ {
+ childAlg->execute();
+ }
+ catch (std::exception&)
+ {
+ g_log.error("Exception thrown while running CropWorkspace as a sub-algorithm");
+ throw;
+ }
- progress(m_fracCompl);
- interruption_point();
+ if ( ! childAlg->isExecuted() )
+ {
+ g_log.error("The CropWorkspace algorithm failed unexpectedly, aborting.");
+ throw std::runtime_error(name() + " failed trying to run CropWorkspace");
+ }
+ return childAlg->getProperty("OutputWorkspace");
+}
+
+/**
+ * Calculate the width of the peak within the given region
+ * @param data_ws The workspace containg the window around the peak
+ * @param prominence The factor that the peak must be above the error to count as a peak
+ * @param peak_x An output vector containing just the X values of the peak data
+ * @param peak_y An output vector containing just the Y values of the peak data
+ * @param peak_e An output vector containing just the E values of the peak data
+ * @returns The width of the peak at half height
+*/
+double GetEi2::calculatePeakWidthAtHalfHeight(MatrixWorkspace_sptr data_ws, const double prominence,
+ std::vector<double> & peak_x, std::vector<double> & peak_y, std::vector<double> & peak_e) const
+{
+ // First create a temporary vector of bin_centre values to work with
+ std::vector<double> Xs(data_ws->readX(0).size());
+ VectorHelper::convertToBinCentre(data_ws->readX(0), Xs);
+ const MantidVec & Ys = data_ws->readY(0);
+ const MantidVec & Es = data_ws->readE(0);
+
+ MantidVec::const_iterator peakIt = std::max_element(Ys.begin(), Ys.end());
+ unsigned int iPeak = peakIt - Ys.begin();
+ double peakY = Ys[iPeak];
+ double peakE = Es[iPeak];
+
+ const unsigned int nxvals = Xs.size();
+
+ //! Find data range that satisfies prominence criterion: im < ipk < ip will be nearest points that satisfy this
+ int im = iPeak-1;
+ for( ; im >= 0; --im )
+ {
+ const double ratio = Ys[im]/peakY;
+ const double ratio_err = std::sqrt( std::pow(Es[im],2) + std::pow(ratio*peakE,2) )/peakY;
+ if ( ratio < (1.0/prominence - m_peak_signif*ratio_err) )
+ {
+ break;
+ }
+ }
+
+ int ip = iPeak+1;
+ for( ; ip < nxvals; ip++ )
+ {
+ const double ratio = Ys[ip]/peakY;
+ const double ratio_err =
+ std::sqrt( std::pow(Es[ip], 2) + std::pow(ratio*peakE, 2) )/peakY;
+ if ( ratio < (1.0/prominence - m_peak_signif*ratio_err) )
+ {
+ break;
}
+ }
+
+ if ( ip == nxvals || im < 0 )
+ {
+ throw std::invalid_argument("No peak found in data that satisfies prominence criterion");
+ }
+
+ // We now have a peak, so can start filling output arguments
+ // Determine extent of peak using derivatives
+ // At this point 1 =< im < ipk < ip =< size(x)
+ // After this section, new values will be given to im, ip that still satisfy these inequalities.
+ //
+ // The algorithm for negative side skipped if im=1; positive side skipped if ip=size(x);
+ // if fails derivative criterion -> ip=size(x) (+ve) im=1 (-ve)
+ // In either case, we deem that the peak has a tail(s) that extend outside the range of x
+
+ if ( ip < nxvals )
+ {
+ double deriv = -1000.0;
+ double error = 0.0;
+ while ( ( ip < nxvals - 1 ) && ( deriv < -m_peak_deriv*error ) )
+ {
+ double dtp = Xs[ip+1] - Xs[ip];
+ double dtm = Xs[ip] - Xs[ip-1];
+ deriv = 0.5*( ((Ys[ip+1] - Ys[ip]) / dtp) + ((Ys[ip] - Ys[ip-1]) / dtm) );
+ error = 0.5*std::sqrt( ( (std::pow(Es[ip+1], 2) + std::pow(Es[ip], 2) ) / std::pow(dtp,2) ) + ((std::pow(Es[ip], 2) + std::pow(Es[ip-1], 2) )/std::pow(dtm,2) )
+ - 2.0*(std::pow(Es[ip], 2) / (dtp*dtm)) );
+ ip++;
+ }
+ ip--;
- /** Implements the Fortran subroute IXFmoments_dataset_2d() from the libISIS
- */
- double GetEi2::getPeakFirstMoments(API::MatrixWorkspace_sptr WS, const double tMin, const double tMax)
+ if (deriv < -error)
{
- // @todo Original FORTRAN uses an unspike algorithm which we don't currently have
- double prominence = 4.0;
- int nvals = WS->dataY(0).size();
- MantidVec centredXs(nvals, 0.0);
+ ip = nxvals -1; // ! derivative criterion not met
+ }
+ }
- for( int i = 1; i <= nvals; i++ )
- {
- centredXs[i - 1] = (WS->readX(0)[i] + WS->readX(0)[i - 1])/2;
- }
+ if (im > 0)
+ {
+ double deriv = 1000.0;
+ double error = 0.0;
+ while ( (im > 0) && (deriv > m_peak_deriv*error) )
+ {
+ double dtp = Xs[im+1] - Xs[im];
+ double dtm = Xs[im] - Xs[im-1];
+ deriv = 0.5*( ((Ys[im+1] - Ys[im]) / dtp) + ( (Ys[im] - Ys[im-1]) / dtm) );
+ error = 0.5*std::sqrt( ( (std::pow(Es[im+1], 2) + std::pow(Es[im], 2) ) / std::pow(dtp, 2) ) + (( std::pow(Es[im], 2) + std::pow(Es[im-1], 2) ) / std::pow(dtm, 2) )
+ - 2.0*std::pow(Es[im], 2)/(dtp*dtm) );
+ im--;
+ }
+ im++;
+ if (deriv > error) im = 0; // derivative criterion not met
+ }
+ double pk_min = Xs[im];
+ double pk_max = Xs[ip];
+ double pk_width = Xs[ip] - Xs[im];
+
+ // Determine background from either side of peak.
+ // At this point, im and ip define the extreme points of the peak
+ // Assume flat background
+ double bkgd = 0.0;
+ double bkgd_range = 0.0;
+ double bkgd_min = std::max(Xs.front(), pk_min - m_bkgd_frac*pk_width);
+ double bkgd_max = std::min(Xs.back(), pk_max + m_bkgd_frac*pk_width);
+
+ if (im > 0)
+ {
+ double bkgd_m, bkgd_err_m;
+ integrate(bkgd_m, bkgd_err_m, Xs, Ys, Es, bkgd_min, pk_min);
+ bkgd = bkgd + bkgd_m;
+ bkgd_range = bkgd_range + (pk_min - bkgd_min);
+ }
- double A_M = 0, c = 0, c_fwhm = 0, w = 0, T_Mean = 0;
- getPeakMean(centredXs, WS->readY(0), WS->readE(0), prominence, A_M, c, c_fwhm, w, T_Mean);
- const double bmin = w/(1.5*NPOINTS);
+ if (ip < nxvals - 1)
+ {
+ double bkgd_p, bkgd_err_p;
+ integrate(bkgd_p, bkgd_err_p, Xs, Ys, Es, pk_max, bkgd_max);
+ bkgd = bkgd + bkgd_p;
+ bkgd_range = bkgd_range + (bkgd_max - pk_max);
+ }
- if (c_fwhm <= 0.0)
+ if ( im > 0 || ip < nxvals - 1 ) //background from at least one side
+ {
+ bkgd = bkgd / bkgd_range;
+ }
+ // Perform moment analysis on the peak after subtracting the background
+ // Fill arrays with peak only:
+ const std::vector<double>::size_type nvalues = ip - im + 1;
+ peak_x.resize(nvalues);
+ std::copy( Xs.begin() + im, Xs.begin() + ip + 1, peak_x.begin());
+ peak_y.resize(nvalues);
+ std::transform( Ys.begin() + im, Ys.begin() + ip + 1, peak_y.begin(), std::bind2nd(std::minus<double>(),bkgd) );
+ peak_e.resize(nvalues);
+ std::copy( Es.begin()+im, Es.begin() + ip + 1, peak_e.begin());
+
+ // FWHH:
+ int ipk_int = iPeak - im; // ! peak position in internal array
+ double hby2 = 0.5*peak_y[ipk_int];
+ int ip1(0), ip2(0);
+ double xp_hh(0);
+
+ const int nyvals = peak_y.size();
+ if (peak_y[nyvals-1] < hby2)
+ {
+ for( int i = ipk_int; i < nyvals; ++i )
+ {
+ if (peak_y[i] < hby2)
{
- throw std::invalid_argument("No peak found, check tMin, tMax and the Monitor index");
+ // after this point the intensity starts to go below half-height
+ ip1 = i-1;
+ break;
}
-
- WS = rebin(WS, tMin, bmin, tMax);
-
- nvals = WS->dataY(0).size();
- centredXs.resize(nvals, 0.0);
-
- for( int i = 1; i <= nvals; i++ )
+ }
+ for ( int i = nyvals-1; i >= ipk_int; --i )
+ {
+ if (peak_y[i]>hby2)
{
- centredXs[i - 1] = (WS->readX(0)[i] + WS->readX(0)[i - 1])/2;
+ ip2 = i+1; // ! point closest to peak after which the intensity is always below half height
+ break;
}
- //! call get moments again, prominence=4 still
- try {
- T_Mean = 0.0;
- getPeakMean(centredXs, WS->readY(0), WS->readE(0), prominence, A_M, c, c_fwhm, w, T_Mean);
- }
- catch(std::invalid_argument &)
+ }
+ xp_hh = peak_x[ip2] + (peak_x[ip1]-peak_x[ip2])*((hby2-peak_y[ip2])/(peak_y[ip1]-peak_y[ip2]));
+ }
+ else
+ {
+ xp_hh = peak_x[nyvals-1];
+ }
+
+ int im1(0), im2(0);
+ double xm_hh(0);
+ if (peak_y[0]<hby2)
+ {
+ for( int i = ipk_int; i >= 0; --i )
+ {
+ if (peak_y[i]<hby2)
{
- prominence = 2.0;
- T_Mean = 0.0;
- getPeakMean(centredXs, WS->readY(0), WS->readE(0), prominence, A_M, c, c_fwhm, w, T_Mean);
+ im1 = i+1; // ! after this point the intensity starts to go below half-height
+ break;
}
-
- if ((c == 0.0) || (w > (0.2*c_fwhm)))
+ }
+ for ( int i=0; i <= ipk_int; ++i )
+ {
+ if (peak_y[i]>hby2)
{
- throw std::invalid_argument("no valid peak found, check initial tMin, tMax and the Monitor index");
+ im2 = i-1; // ! point closest to peak after which the intensity is always below half height
+ break;
}
- return T_Mean;
}
+ xm_hh = peak_x[im2] + (peak_x[im1]-peak_x[im2])*((hby2-peak_y[im2])/(peak_y[im1]-peak_y[im2]));
+ }
+ else
+ {
+ xm_hh = peak_x.front();
+ }
+ return (xp_hh - xm_hh);
+}
+
+/**
+ * Calculate the first moment of the given workspace
+ * @param monitor_ws The workspace containing a single spectrum to work on
+ * @param prominence The factor over the background by which a peak is to be considered a "real" peak
+ */
+double GetEi2::calculateFirstMoment(MatrixWorkspace_sptr monitor_ws, const double prominence)
+{
+ std::vector<double> peak_x, peak_y, peak_e;
+ calculatePeakWidthAtHalfHeight(monitor_ws, prominence, peak_x, peak_y, peak_e);
+
+ // Area
+ double area(0.0), dummy(0.0);
+ double pk_xmin = peak_x.front();
+ double pk_xmax = peak_x.back();
+ integrate(area, dummy, peak_x, peak_y, peak_e, pk_xmin, pk_xmax);
+ // First moment
+ std::transform(peak_y.begin(), peak_y.end(), peak_x.begin(), peak_y.begin(), std::multiplies<double>());
+ double xbar(0.0);
+ integrate(xbar, dummy, peak_x, peak_y, peak_e, pk_xmin, pk_xmax);
+ return xbar / area;
+}
- API::MatrixWorkspace_sptr GetEi2::rebin(API::MatrixWorkspace_sptr WS, const double first, const double width, const double end)
- {
- IAlgorithm_sptr childAlg =
- createSubAlgorithm("Rebin");
+/**
+ * Rebin the workspace using the given bin widths
+ * @param monitor_ws
+ * @param first The minimum value for the new bin range
+ * @param width The new bin width
+ * @param end The maximum value for the new bin range
+ * @returns The rebinned workspace
+*/
+MatrixWorkspace_sptr GetEi2::rebin(MatrixWorkspace_sptr monitor_ws, const double first, const double width, const double end)
+{
+ IAlgorithm_sptr childAlg = createSubAlgorithm("Rebin");
+ childAlg->setProperty("InputWorkspace", monitor_ws);
+ std::ostringstream binParams;
+ binParams << first << "," << width << "," << end;
+ childAlg->setPropertyValue( "Params", binParams.str());
- childAlg->setProperty( "InputWorkspace", WS );
- std::ostringstream binParams;
- binParams << first << "," << width << "," << end;
- childAlg->setPropertyValue( "Params", binParams.str());
+ try
+ {
+ childAlg->execute();
+ }
+ catch (std::exception&)
+ {
+ g_log.error("Exception thrown while running Rebin as a sub-algorithm");
+ throw;
+ }
- try
- {
- childAlg->execute();
- }
- catch (std::exception&)
- {
- g_log.error("Exception thrown while running Regroup as a sub-algorithm");
- throw;
- }
+ if ( ! childAlg->isExecuted() )
+ {
+ g_log.error("The Rebin algorithm failed unexpectedly, aborting.");
+ throw std::runtime_error(name() + " failed trying to run Rebin");
+ }
+ return childAlg->getProperty("OutputWorkspace");
+}
+
+/**
+ * Integrate a point data set
+ * @param integral_val The result of the integral
+ * @param integral_err The error value
+ * @param x The X values
+ * @param s The Y values
+ * @param e The E values
+ * @param xmin The minimum value for the integration
+ * @param xmax The maximum value for the integration
+ */
+void GetEi2::integrate(double & integral_val, double &integral_err, const Mantid::MantidVec &x, const Mantid::MantidVec &s, const Mantid::MantidVec &e,
+ const double xmin, const double xmax) const
+{
+ // MG: Note that this is integration of a point data set from libisis
+ // @todo: Move to Kernel::VectorHelper and improve performance
- if ( ! childAlg->isExecuted() )
- {
- g_log.error("The Rebin algorithm failed unexpectedly, aborting.");
- throw std::runtime_error(name() + " failed trying to run Rebin");
- }
- return childAlg->getProperty("OutputWorkspace");
- }
-
- void GetEi2::getPeakMean(const MantidVec& Xs, const MantidVec& Ys, const MantidVec& Es, const double prominence, double &area, double &c, double &c_fwhm, double &w, double &xbar)
- {
- MantidVec::const_iterator peakIt = std::max_element(Ys.begin(), Ys.end()); //! position of peak
- unsigned int iPeak = peakIt - Ys.begin();
- double peakY = Ys[iPeak];
- double peakE = Es[iPeak];
-
- const unsigned int nxvals = Xs.size();
-
- //! Find data range that satisfies prominence criterion: im < ipk < ip will be nearest points that satisfy this
- int im = iPeak-1;
- for( ; im >= 0; --im )
- {
- const double ratio = Ys[im]/peakY;
- const double ratio_err =
- std::sqrt( std::pow(Es[im],2) + std::pow(ratio*peakE,2) )/peakY;
- if ( ratio < (1.0/prominence - SFAC_PEAK*ratio_err) )
- {
- break;
- }
- }
- if( im < 0 ) im = 0;
-
- int ip = iPeak+1;
- for( ; ip < nxvals; ip++ )
- {
- const double ratio = Ys[ip]/peakY;
- const double ratio_err =
- std::sqrt( std::pow(Es[ip], 2) + std::pow(ratio*peakE, 2) )/peakY;
- if ( ratio < (1.0/prominence - SFAC_PEAK*ratio_err) )
- {
- break;
- }
- }
- if( ip == nxvals ) --ip;
-
- if ( ip < nxvals && im >= 0 ) // ! peak in data
- {
- c = Xs[iPeak];
- }
- else
- {
- throw std::invalid_argument("No peak found in data that satisfies prominence criterion");
- }
- // We now have a peak, so can start filling output arguments
- // Determine extent of peak using derivatives
- // At this point 1 =< im < ipk < ip =< size(x)
- // After this section, new values will be given to im, ip that still satisfy these inequalities.
- //
- // The algorithm for negative side skipped if im=1; positive side skipped if ip=size(x);
- // if fails derivative criterion -> ip=size(x) (+ve) im=1 (-ve)
- // In either case, we deem that the peak has a tail(s) that extend outside the range of x
-
- if ( ip < nxvals )
- {
- double deriv = -1000.0;
- double error = 0.0;
- while ( ( ip < nxvals - 1 ) && ( deriv < -SFAC_DERIV*error ) )
- {
- double dtp = Xs[ip+1] - Xs[ip];
- double dtm = Xs[ip] - Xs[ip-1];
- deriv = 0.5*( ((Ys[ip+1] - Ys[ip]) / dtp) + ((Ys[ip] - Ys[ip-1]) / dtm) );
- error = 0.5*std::sqrt( ( (std::pow(Es[ip+1], 2) + std::pow(Es[ip], 2) ) / std::pow(dtp,2) ) + ((std::pow(Es[ip], 2) + std::pow(Es[ip-1], 2) )/std::pow(dtm,2) )
- - 2.0*(std::pow(Es[ip], 2) / (dtp*dtm)) );
- ip = ip + 1;
- }
- ip = ip - 1;
-
- if (deriv < -error)
- {
- ip = nxvals -1; // ! derivative criterion not met
- }
- }
-
- if (im > 0)
- {
- double deriv = 1000.0;
- double error = 0.0;
- while ( (im > 0) && (deriv > SFAC_DERIV*error) )
- {
- double dtp = Xs[im+1] - Xs[im];
- double dtm = Xs[im] - Xs[im-1];
- deriv = 0.5*( ((Ys[im+1] - Ys[im]) / dtp) + ( (Ys[im] - Ys[im-1]) / dtm) );
- error = 0.5*std::sqrt( ( (std::pow(Es[im+1], 2) + std::pow(Es[im], 2) ) / std::pow(dtp, 2) ) + (( std::pow(Es[im], 2) + std::pow(Es[im-1], 2) ) / std::pow(dtm, 2) )
- - 2.0*std::pow(Es[im], 2)/(dtp*dtm) );
- im = im - 1;
- }
- im = im + 1;
- if (deriv > error) im = 0;// ! derivative criterion not met
- }
- double pk_min = Xs[im];
- double pk_max = Xs[ip];
- double pk_width = Xs[ip] - Xs[im];
-
- // Determine background from either side of peak.
- // At this point, im and ip define the extreme points of the peak
- // Assume flat background
- double bkgd = 0.0;
- double bkgd_range = 0.0;
- double bkgd_min = std::max(Xs.front(), pk_min - BKGD_FAC*pk_width);
- double bkgd_max = std::min(Xs.back(), pk_max + BKGD_FAC*pk_width);
-
- if (im > 0)
- {
- double bkgd_m, bkgd_err_m;
- integrate(bkgd_m, bkgd_err_m, Xs, Ys, Es, bkgd_min, pk_min);
- bkgd = bkgd + bkgd_m;
- bkgd_range = bkgd_range + (pk_min - bkgd_min);
- }
-
- if (ip < nxvals - 1)
- {
- double bkgd_p, bkgd_err_p;
- integrate(bkgd_p, bkgd_err_p, Xs, Ys, Es, pk_max, bkgd_max);
- bkgd = bkgd + bkgd_p;
- bkgd_range = bkgd_range + (bkgd_max - pk_max);
- }
-
- if ( im > 0 || ip < nxvals - 1 ) //background from at least one side
- {
- bkgd = bkgd / bkgd_range;
- }
- // Perform moment analysis on the peak after subtracting the background
- // Fill arrays with peak only:
- std::vector<double> xint(ip-im+1);
- std::copy( Xs.begin()+im, Xs.begin()+ip+1, xint.begin());
- std::vector<double> yint(ip-im+1);
- std::transform( Ys.begin()+im, Ys.begin()+ip+1, yint.begin(), std::bind2nd(std::minus<double>(),bkgd) );
- std::vector<double> eint(ip-im+1);
- std::copy( Es.begin()+im, Es.begin()+ip+1, eint.begin());
-
- // ! FWHH:
- int ipk_int = iPeak-im; // ! peak position in internal array
- double hby2 = 0.5*yint[ipk_int];
- int ip1(0), ip2(0);
- double xp_hh(0);
-
- const int nyvals = yint.size();
- if (yint[nyvals-1]<hby2)
- {
- for( int i = ipk_int; i < nyvals; ++i )
- {
- if (yint[i] < hby2)
- {
- // after this point the intensity starts to go below half-height
- ip1 = i-1;
- break;
- }
- }
- for ( int i = nyvals-1; i >= ipk_int; --i )
- {
- if (yint[i]>hby2)
- {
- ip2 = i+1; // ! point closest to peak after which the intensity is always below half height
- break;
- }
- }
- xp_hh = xint[ip2] + (xint[ip1]-xint[ip2])*((hby2-yint[ip2])/(yint[ip1]-yint[ip2]));
- }
- else
- {
- xp_hh = xint[nyvals-1];
- }
+ MantidVec::const_iterator lowit = std::lower_bound(x.begin(), x.end(), xmin);
+ MantidVec::difference_type ml = std::distance(x.begin(),lowit);
+ MantidVec::const_iterator highit = std::upper_bound(lowit,x.end(), xmax);
+ MantidVec::difference_type mu = std::distance(x.begin(),highit);
+ if( mu > 0 ) --mu;
+
+ MantidVec::size_type nx(x.size());
+ if( mu < ml )
+ {
+ //special case of no data points in the integration range
+ unsigned int ilo = std::max<unsigned int>(ml-1,0);
+ unsigned int ihi = std::min<unsigned int>(mu+1,nx);
+ double fraction = (xmax - xmin)/(x[ihi] - x[ilo]);
+ integral_val = 0.5 * fraction *
+ ( s[ihi]*((xmax - x[ilo]) + (xmin - x[ilo])) + s[ilo]*((x[ihi] - xmax)+(x[ihi] - xmin)) );
+ double err_hi = e[ihi]*((xmax - x[ilo]) + (xmin - x[ilo]));
+ double err_lo = e[ilo]*((x[ihi] - xmax)+(x[ihi] - xmin));
+ integral_err = 0.5*fraction*std::sqrt( err_hi*err_hi + err_lo*err_lo );
+ return;
+ }
+
+ double x1eff(0.0), s1eff(0.0), e1eff(0.0);
+ if( ml > 0 )
+ {
+ x1eff = (xmin*(xmin - x[ml-1]) + x[ml-1]*(x[ml] - xmin))/(x[ml] - x[ml-1]);
+ double fraction = (x[ml]-xmin) / ((x[ml]-x[ml-1]) + (xmin - x[ml-1]));
+ s1eff = s[ml-1]* fraction;
+ e1eff = e[ml-1]* fraction;
+ }
+ else
+ {
+ x1eff = x[ml];
+ s1eff = 0.0;
+ }
- int im1(0), im2(0);
- double xm_hh(0);
- if (yint[0]<hby2)
- {
- for ( int i=ipk_int; i >= 0; --i )
- {
- if (yint[i]<hby2)
- {
- im1 = i+1; // ! after this point the intensity starts to go below half-height
- break;
- }
- }
- for ( int i=0; i <= ipk_int; ++i )
- {
- if (yint[i]>hby2)
- {
- im2 = i-1; // ! point closest to peak after which the intensity is always below half height
- break;
- }
- }
- xm_hh = xint[im2] + (xint[im1]-xint[im2])*((hby2-yint[im2])/(yint[im1]-yint[im2]));
- }
- else
+ double xneff(0.0), sneff(0.0), eneff;
+ if( mu < nx - 1)
+ {
+ xneff = (xmax*(x[mu+1]-xmax) + x[mu+1]*(xmax - x[mu]))/(x[mu+1] - x[mu]);
+ const double fraction = (xmax - x[mu])/((x[mu+1] - x[mu]) + (x[mu+1]-xmax));
+ sneff = s[mu+1]*fraction;
+ eneff = e[mu+1]*fraction;
+ }
+ else
+ {
+ xneff = x.back();
+ sneff = 0.0;
+ }
+
+ //xmin -> x[ml]
+ integral_val = (x[ml] - x1eff)*(s[ml] + s1eff);
+ integral_err = e1eff*(x[ml] - x1eff);
+ integral_err *= integral_err;
+
+ if( mu == ml )
+ {
+ double ierr = e[ml]*(xneff - x1eff);
+ integral_err += ierr * ierr;
+ }
+ else if( mu == ml + 1 )
+ {
+ integral_val += (s[mu] + s[ml])*(x[mu] - x[ml]);
+ double err_lo = e[ml]*(x[ml+1] - x1eff);
+ double err_hi = e[mu]*(x[mu-1] - xneff);
+ integral_err += err_lo*err_lo + err_hi*err_hi;
+ }
+ else
+ {
+ for( int i = ml; i < mu; ++i )
+ {
+ integral_val += (s[i+1] + s[i])*(x[i+1] - x[i]);
+ if( i < mu - 1 )
{
- xm_hh = xint[0];
+ double ierr = e[i+1]*(x[i+2] - x[i]);
+ integral_err += ierr * ierr;
}
-
- c_fwhm = 0.5*(xp_hh + xm_hh);
- w = xp_hh - xm_hh;
-
- // ! area:
- double dummy;
- integrate(area, dummy, xint, yint, eint, pk_min, pk_max);
- // ! first moment:
- std::transform(yint.begin(), yint.end(), xint.begin(), yint.begin(), std::multiplies<double>());
- integrate(xbar, dummy, xint, yint, eint, pk_min, pk_max);
- xbar = xbar / area;
-
- // look out for user cancel messgages
- advanceProgress(FIT_PEAK);
}
-
- void GetEi2::integrate(double & integral_val, double &integral_err, const MantidVec &x, const MantidVec &s, const MantidVec &e, const double xmin, const double xmax)
- {
- // MG: Note that this is integration of a point data set from libisis
- // @todo: Move to Kernel::VectorHelper and improve performance
+ }
+
+ //x[mu] -> xmax
+ integral_val += (xneff - x[mu])*(s[mu] + sneff);
+ double err_tmp = eneff*(xneff - x[mu]);
+ integral_err += err_tmp*err_tmp;
+
+ integral_val *= 0.5;
+ integral_err = 0.5*sqrt(integral_err);
+}
+
+/**
+ * Offset the bins on the input workspace by the calculated time of the first monitor peak
+ */
+void GetEi2::applyBinOffset()
+{
+ IAlgorithm_sptr child_alg = createSubAlgorithm("ChangeBinOffset");
+ child_alg->setProperty("InputWorkspace", m_input_ws);
+ child_alg->setProperty("Offset", -m_peak1_pos.second);
- MantidVec::const_iterator lowit = std::lower_bound(x.begin(), x.end(), xmin);
- MantidVec::difference_type ml = std::distance(x.begin(),lowit);
- MantidVec::const_iterator highit = std::upper_bound(lowit,x.end(), xmax);
- MantidVec::difference_type mu = std::distance(x.begin(),highit);
- if( mu > 0 ) --mu;
-
- MantidVec::size_type nx(x.size());
- if( mu < ml )
- {
- //special case of no data points in the integration range
- unsigned int ilo = std::max<unsigned int>(ml-1,0);
- unsigned int ihi = std::min<unsigned int>(mu+1,nx);
- double fraction = (xmax - xmin)/(x[ihi] - x[ilo]);
- integral_val = 0.5 * fraction *
- ( s[ihi]*((xmax - x[ilo]) + (xmin - x[ilo])) + s[ilo]*((x[ihi] - xmax)+(x[ihi] - xmin)) );
- double err_hi = e[ihi]*((xmax - x[ilo]) + (xmin - x[ilo]));
- double err_lo = e[ilo]*((x[ihi] - xmax)+(x[ihi] - xmin));
- integral_err = 0.5*fraction*std::sqrt( err_hi*err_hi + err_lo*err_lo );
- return;
- }
-
- double x1eff(0.0), s1eff(0.0), e1eff(0.0);
- if( ml > 0 )
- {
- x1eff = (xmin*(xmin - x[ml-1]) + x[ml-1]*(x[ml] - xmin))/(x[ml] - x[ml-1]);
- double fraction = (x[ml]-xmin) / ((x[ml]-x[ml-1]) + (xmin - x[ml-1]));
- s1eff = s[ml-1]* fraction;
- e1eff = e[ml-1]* fraction;
- }
- else
- {
- x1eff = x[ml];
- s1eff = 0.0;
- }
-
- double xneff(0.0), sneff(0.0), eneff;
- if( mu < nx - 1)
- {
- xneff = (xmax*(x[mu+1]-xmax) + x[mu+1]*(xmax - x[mu]))/(x[mu+1] - x[mu]);
- const double fraction = (xmax - x[mu])/((x[mu+1] - x[mu]) + (x[mu+1]-xmax));
- sneff = s[mu+1]*fraction;
- eneff = e[mu+1]*fraction;
- }
- else
- {
- xneff = x.back();
- sneff = 0.0;
- }
-
- //xmin -> x[ml]
- integral_val = (x[ml] - x1eff)*(s[ml] + s1eff);
- integral_err = e1eff*(x[ml] - x1eff);
- integral_err *= integral_err;
+ try
+ {
+ child_alg->execute();
+ }
+ catch(std::exception& e)
+ {
+ g_log.error() << e.what() << "\n";
+ g_log.error() << "Error executing ChangeBinOffset algorithm. Input bins kept.\n";
+ }
+ m_input_ws = child_alg->getProperty("OutputWorkspace");
+}
- if( mu == ml )
- {
- double ierr = e[ml]*(xneff - x1eff);
- integral_err += ierr * ierr;
- }
- else if( mu == ml + 1 )
- {
- integral_val += (s[mu] + s[ml])*(x[mu] - x[ml]);
- double err_lo = e[ml]*(x[ml+1] - x1eff);
- double err_hi = e[mu]*(x[mu-1] - xneff);
- integral_err += err_lo*err_lo + err_hi*err_hi;
- }
- else
- {
- for( int i = ml; i < mu; ++i )
- {
- integral_val += (s[i+1] + s[i])*(x[i+1] - x[i]);
- if( i < mu - 1 )
- {
- double ierr = e[i+1]*(x[i+2] - x[i]);
- integral_err += ierr * ierr;
- }
- }
- }
-
- //x[mu] -> xmax
- integral_val += (xneff - x[mu])*(s[mu] + sneff);
- double err_tmp = eneff*(xneff - x[mu]);
- integral_err += err_tmp*err_tmp;
-
- integral_val *= 0.5;
- integral_err = 0.5*sqrt(integral_err);
- }
-
- /** Get the kinetic energy of a neuton in joules given it speed using E=mv^2/2
- * @param speed the instantanious speed of a neutron in metres per second
- * @return the energy in joules
- */
- double GetEi2::neutron_E_At(double speed) const
- {
- // E_KE = mv^2/2
- return PhysicalConstants::NeutronMass*speed*speed/(2);
- }
-
- /// Update the percentage complete estimate assuming that the algorithm has completed a task with estimated RunTime toAdd
- void GetEi2::advanceProgress(double toAdd)
- {
- m_fracCompl += toAdd;
- progress(m_fracCompl);
- // look out for user cancel messgages
- interruption_point();
- }
-
-
- } // namespace Algorithms
-} // namespace Mantid
+/**
+ * Move the source of neutrons to the position defined by the first input monitor
+ */
+void GetEi2::moveNeutronSource()
+{
+ Geometry::IDetector_sptr mon_1;
+ try
+ {
+ mon_1 = m_input_ws->getDetector(m_peak1_pos.first);
+ }
+ catch(std::exception&)
+ {
+ g_log.warning() << "Error retrieving detector corresponding workspace index " << m_peak1_pos.first << ". Cannot move component.";
+ return;
+ }
+
+ Geometry::IObjComponent_sptr neutron_src = m_input_ws->getInstrument()->getSource();
+ const Geometry::V3D mon1_pos = mon_1->getPos();
+
+ IAlgorithm_sptr child_alg = createSubAlgorithm("MoveInstrumentComponent");
+ child_alg->setProperty("Workspace", m_input_ws);
+ child_alg->setProperty("ComponentName", neutron_src->getName());
+ child_alg->setProperty("X", mon1_pos.X());
+ child_alg->setProperty("Y", mon1_pos.Y());
+ child_alg->setProperty("Z", mon1_pos.Z());
+ child_alg->setProperty("RelativePosition", "0");
+
+ try
+ {
+ child_alg->execute();
+ }
+ catch(std::exception& e)
+ {
+ g_log.error() << e.what() << "\n";
+ g_log.error() << "Error executing MoveInstrumentComponent as sub algorithm. Neutron source unmoved\n";
+ }
+ m_input_ws = child_alg->getProperty("Workspace");
+
+}
+
+/**
+ * Store the value of Ei wihin the log data on the Sample object
+ * @param ei The value of the incident energy of the neutron
+ */
+void GetEi2::storeEi(const double ei) const
+{
+ Property *incident_energy = new PropertyWithValue<double>("Ei",ei,Direction::Input);
+ m_input_ws->mutableSample().addLogData(incident_energy);
+}
\ No newline at end of file
diff --git a/Code/Mantid/DataHandling/inc/MantidDataHandling/LoadDetectorInfo.h b/Code/Mantid/DataHandling/inc/MantidDataHandling/LoadDetectorInfo.h
index 7bcc9f2c5adb8628dded9df33e768c7483572851..2b6612df728657ae1f39693354477782d531eefb 100644
--- a/Code/Mantid/DataHandling/inc/MantidDataHandling/LoadDetectorInfo.h
+++ b/Code/Mantid/DataHandling/inc/MantidDataHandling/LoadDetectorInfo.h
@@ -74,10 +74,6 @@ private:
};
/// will store a pointer to the user selected workspace
DataObjects::Workspace2D_sptr m_workspace;
- /// the instrument with in the user selected workspace
- API::IInstrument_sptr m_instrument;
- /// the map that stores additional properties for detectors
- Geometry::ParameterMap *m_paraMap;
/// number of histograms in the workspace
int m_numHists;
/// the detector IDs that are monitors, according to the raw file
diff --git a/Code/Mantid/DataHandling/src/LoadDetectorInfo.cpp b/Code/Mantid/DataHandling/src/LoadDetectorInfo.cpp
index 4178d4ca48f751e147e74fad784273b0b32e4f58..81fab673388342516d02db4716dfac3b09b8582e 100644
--- a/Code/Mantid/DataHandling/src/LoadDetectorInfo.cpp
+++ b/Code/Mantid/DataHandling/src/LoadDetectorInfo.cpp
@@ -27,7 +27,7 @@ const LoadDetectorInfo::detectDatForm
LoadDetectorInfo::MAPS_MER_TYPE(14, 11, 12);
/// Empty default constructor
LoadDetectorInfo::LoadDetectorInfo() : Algorithm(),
- m_workspace(), m_instrument(), m_paraMap(NULL), m_numHists(-1), m_monitors(),
+ m_workspace(), m_numHists(-1), m_monitors(),
m_monitorXs(), m_commonXs(false), m_monitOffset(UNSETOFFSET), m_error(false),
m_FracCompl(0.0)
{
@@ -64,8 +64,6 @@ void LoadDetectorInfo::exec()
{
// get the infomation that will be need from the user selected workspace, assume it exsists because of the validator in init()
m_workspace = getProperty("Workspace");
- m_instrument = m_workspace->getInstrument();
- m_paraMap = &(m_workspace->instrumentParameters());
m_numHists = m_workspace->getNumberHistograms();
// when we change the X-values will care take to maintain sharing. I have only implemented maintaining sharing where _all_ the arrays are initiall common
m_commonXs = WorkspaceHelpers::sharedXData(m_workspace);
@@ -382,31 +380,43 @@ void LoadDetectorInfo::readRAW(const std::string& fName)
*/
void LoadDetectorInfo::setDetectorParams(const detectorInfo ¶ms, detectorInfo &change)
{
- IComponent* baseComp = NULL;
+ Geometry::IDetector_sptr det;
try
- {// there are often exceptions caused by detectors not being listed in the instrument file
- Geometry::IDetector_sptr det = m_instrument->getDetector(params.detID);
- baseComp = det->getComponent();
+ {
+ det = m_workspace->getInstrument()->getDetector(params.detID);
}
catch( std::runtime_error &e)
- {// when detector information is incomplete there are a range of exceptions you can get, I think none of which should kill the algorithm
+ {
throw Exception::NotFoundError(e.what(), params.detID);
}
+ Geometry::ParameterMap &pmap = m_workspace->instrumentParameters();
+ IComponent* comp = det->getComponent();
// set the detectors pressure, first check if it already has a setting, if not add it
- Parameter_sptr setting = m_paraMap->get(baseComp, "3He(atm)");
- if (setting) setting->set(params.pressure);
+ Parameter_sptr setting = pmap.get(comp, "3He(atm)");
+ if (setting)
+ {
+ setting->set(params.pressure);
+ }
else
- m_paraMap->add("double", baseComp, "3He(atm)", params.pressure);
+ {
+ pmap.add("double", comp, "3He(atm)", params.pressure);
+ }
- setting = m_paraMap->get(baseComp, "wallT(m)");
- if (setting) setting->set(params.wallThick);
+ setting = pmap.get(comp, "wallT(m)");
+ if (setting)
+ {
+ setting->set(params.wallThick);
+ }
else
- m_paraMap->add("double", baseComp, "wallT(m)", params.wallThick);
+ {
+ pmap.add("double", comp, "wallT(m)", params.wallThick);
+ }
// this operation has been successful if we are here, the following infomation is usefull for logging
change = params;
}
+
/** Decides if the bin boundaries for all non-monitor spectra will be the same and runs the appropiate
* function. It is possible for this function to set all detectors spectra X-values when they shouldn't
* @param lastOffset the delay time of the last detector is only used if differentDelays is false or if detectID and delays are leave empty e.g. when we use common bins
diff --git a/Code/Mantid/Kernel/inc/MantidKernel/VectorHelper.h b/Code/Mantid/Kernel/inc/MantidKernel/VectorHelper.h
index aeef4e2e100826d5d28931f6556c112756c9fd62..324b7cfb6e134ea42726c4b2eebf200669d1c9f6 100644
--- a/Code/Mantid/Kernel/inc/MantidKernel/VectorHelper.h
+++ b/Code/Mantid/Kernel/inc/MantidKernel/VectorHelper.h
@@ -51,6 +51,9 @@ namespace VectorHelper
void DLLExport rebinHistogram(const std::vector<double>& xold, const std::vector<double>& yold, const std::vector<double>& eold,
const std::vector<double>& xnew, std::vector<double>& ynew, std::vector<double>& enew,bool addition);
+ /// Convert an array of bin boundaries to bin centre values.
+ void DLLExport convertToBinCentre(const std::vector<double> & bin_edges, std::vector<double> & bin_centres);
+
//! Functor used for computing the sum of the square values of a vector, using the accumulate algorithm
template <class T> struct SumGaussError: public std::binary_function<T,T,T>
{
@@ -118,6 +121,17 @@ namespace VectorHelper
}
};
+ /// A binary functor to compute the simple average of 2 numbers
+ template <class T> struct SimpleAverage: public std::binary_function<T,T,T>
+ {
+ SimpleAverage() {}
+ /// Return the average of the two arguments
+ T operator()(const T & x, const T & y) const
+ {
+ return 0.5 * (x + y);
+ }
+ };
+
} // namespace VectorHelper
} // namespace Kernel
} // namespace Mantid
diff --git a/Code/Mantid/Kernel/src/VectorHelper.cpp b/Code/Mantid/Kernel/src/VectorHelper.cpp
index 7f1f96fe3f4c55dab5e3308d12811a4ab27ffc36..8459d8fdbd59fee3066a5b0b54c3fcbf6f2ffbeb 100644
--- a/Code/Mantid/Kernel/src/VectorHelper.cpp
+++ b/Code/Mantid/Kernel/src/VectorHelper.cpp
@@ -3,6 +3,7 @@
#include "MantidKernel/VectorHelper.h"
#include <algorithm>
+#include <numeric>
#include <iostream>
namespace Mantid
@@ -283,6 +284,27 @@ void rebinHistogram(const std::vector<double>& xold, const std::vector<double>&
return;
}
+/**
+ * Convert the given set of bin boundaries into bin centre values
+ * @param bin_edges A vector of values specifying bin boundaries
+ * @param bin_centres An output vector of bin centre values.
+*/
+void convertToBinCentre(const std::vector<double> & bin_edges, std::vector<double> & bin_centres)
+{
+ const std::vector<double>::size_type npoints = bin_edges.size();
+ if( bin_centres.size() != npoints )
+ {
+ bin_centres.resize(npoints);
+ }
+
+ // The custom binary function modifies the behaviour of the algorithm to compute the average of
+ // two adjacent bin boundaries
+ std::adjacent_difference(bin_edges.begin(), bin_edges.end(), bin_centres.begin(), SimpleAverage<double>());
+ // The algorithm copies the first element of the input to the first element of the output so we need to
+ // remove the first element of the output
+ bin_centres.erase(bin_centres.begin());
+}
+
} // End namespace VectorHelper
} // End namespace Kernel
} // End namespace Mantid
diff --git a/Code/Mantid/PythonAPI/scripts/Excitations/CommonFunctions.py b/Code/Mantid/PythonAPI/scripts/Excitations/CommonFunctions.py
index 933f5b1106461a6b790e99fcb6665b449e1fab30..7f1377739412ede663278959a78f04788d85c7c7 100644
--- a/Code/Mantid/PythonAPI/scripts/Excitations/CommonFunctions.py
+++ b/Code/Mantid/PythonAPI/scripts/Excitations/CommonFunctions.py
@@ -31,41 +31,45 @@ def stringToList(commaSeparated):
return theList
#sum all the workspaces, when the workspaces are not summed single input files are specified in this file and the final Python script is made of many copies of this file
-def sumWorkspaces(total, instrument, runNumbers):
+def sumWorkspaces(total, prefix, runNumbers):
if len(runNumbers) > 1:
tempWS = 'CommonFuncs_sumWorkspaces_tempory'
- for toAdd in runNumbers[ 1 : ] :
- instrument.loadRun(toAdd, tempWS) #workspaces are overwriten to save memory
+ for toAdd in runNumbers[1:]:
+ loadRun(prefix, toAdd, tempWS)
Plus(total, tempWS, total)
- mantid.deleteWorkspace(tempWS) #from www.mantidproject.org/MantidPyFramework
+ mantid.deleteWorkspace(tempWS)
#-- Functions to do with input files
# uses the extension to decide whether use LoadNexus or LoadRaw
def LoadNexRaw(filename, workspace):
filename = filename.strip() #remove any white space from the front or end of the name
# this removes everything after the last '.' (rpartition always returns three strings)
- extension = filename.rpartition('.')[2]
- if (extension.lower() == 'nxs') :
+ partitions = filename.split('.')
+ if len(partitions) != 2:
+ raise RuntimeError('Incorrect filename format encountered "' + filename + '"')
+ extension = filename.split('.')[1]
+ if (extension.lower() == 'raw'):
+ loader = LoadRaw(filename, workspace)
+ elif (extension.lower() == 'nxs'):
#return the first property from the algorithm, which for LoadNexus is the output workspace
- return LoadNexus(filename, workspace)[0]
- if (extension.lower() == 'raw') :
- return LoadRaw(filename, workspace)[0]
- #we shouldn't get to here, the function should have returned by now
- raise Exception("Could not find a load function for file "+filename+", *.raw and *.nxs accepted")
+ loader = LoadNexus(filename, workspace)
+ else:
+ raise Exception("Could not find a load function for file "+filename+", *.raw and *.nxs accepted")
+ return loader.workspace(),loader.getPropertyValue("Filename")
# guess the filename from run number using the instrument code
def getFileName(instrumentPref, runNumber):
runNumber = str(runNumber)
- try :
- number = int(runNumber, 10)
- except TypeError :
+ try:
+ number = int(runNumber)
+ except ValueError:
# means we weren't passed a number assume it is a valid file name and return it unprocessed
return runNumber
#only raw files are supported at the moment
return instrumentPref + runNumber + '.raw'
-def loadRun(instru, runNum, workspace):
- return LoadNexRaw(getFileName(instru, runNum), workspace)
+def loadRun(prefix, runNum, workspace):
+ return LoadNexRaw(getFileName(prefix, runNum), workspace)
def loadMask(MaskFilename):
inFile = open(MaskFilename)
@@ -75,7 +79,6 @@ def loadMask(MaskFilename):
if len(numbers) > 0 : # ignore empty lines
if numbers[0][0].isdigit() : # any non-numeric character at the start of the line marks a comment, check the first character of the first word
for specNumber in numbers :
-
if specNumber == '-' :
spectraList[len(spectraList)-1] = spectraList + '-' #if there is a hyphen we don't need commas
else : spectraList = spectraList + "," + specNumber
@@ -84,7 +87,7 @@ def loadMask(MaskFilename):
mantid.sendLogMessage('Only comment lines found in mask file '+MaskFilename)
return ''
return spectraList[1:] #return everything after the very first comma we added in the line above
-
+
def getRunName(path):
# get the string after the last /
filename = path.split('/')
@@ -95,10 +98,14 @@ def getRunName(path):
# remove the last '.' and everything after it i.e. the extension. If there is not extension this just returns the whole thing
return filename.rpartition('.')[0]
+def loadRun(prefix, runNum, workspace):
+ return LoadNexRaw(getFileName(prefix,runNum), workspace)
+
#-- Holds data about the defaults used for diferent instruments (MARI, MAPS ...)
class defaults:
# set the defaults for a default machine. These default values for defaults won't work and so they must be overriden by the correct values for the machine when this is run
- def __init__(self, background_range=(-1.0,-1.0), normalization='not set', instrument_pref='', white_beam_integr=(-1.0,-1.0), scale_factor=1, monitor1_integr=(-1.0e5, -1.0e5), white_beam_scale = 1.0):
+ def __init__(self, background_range=(-1.0,-1.0), normalization='not set', instrument_pref='', white_beam_integr=(-1.0,-1.0), \
+ scale_factor=1, monitor1_integr=(-1.0e5, -1.0e5), white_beam_scale = 1.0, getei_monitors = (-1,-1)):
self.background_range=background_range
self.normalization=normalization
self.instrument_pref=instrument_pref
@@ -106,6 +113,7 @@ class defaults:
self.scale_factor=scale_factor
self.monitor1_integr=monitor1_integr
self.white_beam_scale = white_beam_scale
+ self.getei_monitors = getei_monitors
# guess the filename from run number assuming global_getFileName_instrument_pref is setup
def getFileName(self, runNumber):
@@ -116,5 +124,3 @@ class defaults:
return runNumber
return getFileName(self.instrument_pref, number)
- def loadRun(self, runNum, workspace):
- return LoadNexRaw(self.getFileName(runNum), workspace)
diff --git a/Code/Mantid/PythonAPI/scripts/Excitations/ConversionLib.py b/Code/Mantid/PythonAPI/scripts/Excitations/ConversionLib.py
index eef65b2d8e9ce6401a2cabe4ba42f718c8f5fddf..934d9f73edcf484e2f7b027379fc01aae78cd8fc 100644
--- a/Code/Mantid/PythonAPI/scripts/Excitations/ConversionLib.py
+++ b/Code/Mantid/PythonAPI/scripts/Excitations/ConversionLib.py
@@ -1,165 +1,134 @@
+"""Excitations analysis module
+"""
from mantidsimple import *
import CommonFunctions as common
-# these are the defaults for different instruments
-import ExcitDefaults
import math
-def NormaliseTo(reference, WS, fromTOF, instrument=common.defaults()):
- if (reference == 'monitor-monitor 1') :
- if instrument.monitor1_integr[0] <= -1e5 or instrument.monitor1_integr[1] <= -1e5 :
- raise Exception('Can\'t normalize to monitor1 without instrument defaults information')
- min = instrument.monitor1_integr[0] - fromTOF
- max = instrument.monitor1_integr[1] - fromTOF
- NormaliseToMonitor( InputWorkspace=WS, OutputWorkspace=WS, MonitorSpectrum=1, IntegrationRangeMin=min, IntegrationRangeMax=max,IncludePartialBins=True)
-
- elif reference == 'monitor-monitor 2' :
-
- PEAKWIDTH = 4.0/100
- # the origin of time was set to the location of the peak so to get it's area integrate around the origin, PEAKWIDTH fraction of total time of flight
- min = (-PEAKWIDTH/2)*20000
- max = (PEAKWIDTH/2)*20000
-
- NormaliseToMonitor( InputWorkspace=WS, OutputWorkspace=WS, MonitorSpectrum=2, IntegrationRangeMin=min, IntegrationRangeMax=max)
+# these are the defaults for different instruments. Need to be moved to parameter file
+import ExcitDefaults
- elif reference == 'protons (uAh)' :
- NormaliseByCurrent( InputWorkspace=WS, OutputWorkspace=WS )
+def mono_sample(inst_prefix, run_nums, Ei, d_rebin, wbrf, getEi=True, back='', norma='', det_map='', det_mask='', output_name = 'mono_sample_temporyWS') :
+ """
+ Calculate Ei and detector offsets for a mono-chromatic run.
+ """
+ run_nums = common.listToString(run_nums).split(',')
+ result_ws, det_eff_file = common.loadRun(inst_prefix,run_nums[0], output_name)
+ if result_ws.isGroup():
+ raise RunTimeError("Workspace groups are not supported here")
+
+ instrument = result_ws.getInstrument()
+ # The final workspace will take the X-values and instrument from the first workspace and so we don't need to rerun ChangeBinOffset(), MoveInstrumentComponent(), etc.
+ common.sumWorkspaces(result_ws, inst_prefix, run_nums)
+ if det_eff_file != '':
+ LoadDetectorInfo(result_ws, det_eff_file)
+
+ Ei, mon1_peak = calculateEi(result_ws, Ei, getEi,instrument)
+ bin_offset = -mon1_peak
+
+ if back != 'noback':
+ TOFLow = instrument.getNumberParameter("bkgd-range-min")[0]
+ TOFHigh = instrument.getNumberParameter("bkgd-range-max")[0]
+ # Remove the count rate seen in the regions of the histograms defined as the background regions, if the user defined a region
+ ConvertToDistribution(result_ws)
+ FlatBackground(result_ws, result_ws, TOFLow + bin_offset, TOFHigh + bin_offset, '', 'Mean')
+ ConvertFromDistribution(result_ws)
- elif reference != 'no normalization' :
- raise Exception('Normalisation scheme ' + reference + ' not found. It must be one of monitor, current, peak or none')
+ # deals with normalize to monitor (e.g. norm = 'monitor-monitor 1'), current (if norm = 'protons (uAh)'), etc.
+ NormaliseTo(norma, result_ws, bin_offset, instrument)
-def NormaliseToWhiteBeam(instr, WBRun, toNorm, mapFile, detMask, prevNorm) :
- theNorm = "_ETrans_norm_tempory_WS"
- try:
- pNorm = common.LoadNexRaw(instr.getFileName(WBRun), theNorm)
- NormaliseTo(prevNorm, pNorm, 0, instr)
+ ConvertUnits(result_ws, result_ws, 'DeltaE', 'Direct', Ei, AlignBins=0)
- ConvertUnits(pNorm, pNorm, "Energy", AlignBins=0)
-
- #this both integrates the workspace into one bin spectra and sets up common bin boundaries for all spectra
- low = instr.white_beam_integr[0]
- hi = instr.white_beam_integr[1]
- delta = 2.0*(hi - low)
- Rebin(pNorm, pNorm, str(low)+','+str(delta)+','+str(hi) )
-
- if detMask != '' :
- MaskDetectors(pNorm, SpectraList=detMask)
-
- if mapFile!= "" :
- GroupDetectors(pNorm, pNorm, mapFile, KeepUngroupedSpectra=0)
-
- toNorm /= pNorm
-
- finally:
- mantid.deleteWorkspace(theNorm)
-
-# returns the background range from the string range, which could be two numbers separated by a comma or could be empty and then the default values in instrument are used. Another possiblity is that background correction isn't going to be done, then we don't have to do anything
-# throws if the string range is in the wrong format
-def getBackRange(range, instrument) :
- if range == 'noback' : return ['noback', 'noback']
- # load the default time of flight values that delimit the background region
- TOFLow = instrument.background_range[0]
- TOFHigh = instrument.background_range[1]
- if range != '' :
- range = common.listToString(range)
- strings = range.split(',')
- TOFLow = float(strings[0])
- TOFHigh = float(strings[1])
- return [TOFLow, TOFHigh]
-
-## return specified normalisation method or the default for the instrument if '' was passed
-def getNorm(instrument, normMethod):
- if normMethod == '' : return instrument.normalization
- return normMethod
-
-###--Read in user set parameters or defaults values---------------
-def retrieveSettings(instrum, back, runs):
- instru = ExcitDefaults.loadDefaults(instrum)
+ if d_rebin != '':
+ Rebin(result_ws, result_ws, common.listToString(d_rebin))
+
+ if det_eff_file != '':
+ DetectorEfficiencyCor(result_ws, result_ws, Ei)
+
+ if det_mask != '' :
+ MaskDetectors(Workspace=result_ws, SpectraList=common.listToString(det_mask))
- run_nums = common.listToString(runs).split(',')
- try:
- [TOFLow, TOFHigh] = getBackRange(back, instru)
- except:
- raise Exception('The background range must be specified as a string\n of two numbers separated by a comma specifing the range\nor noback to disable background correction')
+ if det_map != '':
+ GroupDetectors(result_ws, result_ws, det_map, KeepUngroupedSpectra=0)
+ ConvertToDistribution(result_ws)
- return (instru, TOFLow, TOFHigh, run_nums)
-
-def findPeaksOffSetTimeAndEi(workS, E_guess, getEi=True, detInfoFile=''):
- if detInfoFile != '' :
- LoadDetectorInfo(workS, detInfoFile) #for details see www.mantidproject.org/LoadDetectorInfo
-
- runGetEi = str(getEi).lower() != 'fixei' and str(getEi).lower() != 'false'
- if not runGetEi :
- mtd.sendLogMessage('Getting peak times (the GetEi() energy will be overriden by '+str(E_guess)+' meV)')
-
- # Get incident energy
- GetEiData = GetEi(workS, 2, 3, E_guess)
+ NormaliseToWhiteBeam(wbrf, result_ws, det_map, det_mask, norma, inst_prefix, instrument)
- if runGetEi : Ei = GetEiData.getPropertyValue('IncidentEnergy')
- else : Ei = E_guess
+ # Overall scale factor
+ scale_factor = instrument.getNumberParameter("scale-factor")[0]
+ result_ws *= scale_factor
- offSetTime=float(GetEiData.getPropertyValue('FirstMonitorPeak')) # set the origin of time to be when the neutrons arrive at the first GetEi monitor (often called monitor 2)
- ChangeBinOffset(workS, workS, -offSetTime)
-
- #??STEVES?? MARI specfic code GetEiData.FirstMonitor ?
- p0=workS.getDetector(1).getPos() # set the source to be where the neutrons were at the origin of time (that first GetEi monitor)
- mtd.sendLogMessage('Adjusting the X-values so that the zero of time is when the peak was registered at the new origin , which monitor 2')
- src = workS.getInstrument().getSource().getName() #this name could be 'Moderator', 'undulator', etc.
- MoveInstrumentComponent(workS, src, X=p0.getX() , Y=p0.getY() ,Z=p0.getZ(), RelativePosition=False)
- return Ei, offSetTime
+ return result_ws
+
+def calculateEi(input_ws, E_guess, getEi, instrument):
+ """
+ Calculate incident energy of neutrons
+ """
+ getEi = str(getEi).lower()
+ if getEi == 'fixei' or getEi == 'false':
+ fixei = True
+ else:
+ fixei = False
-# a workspace name that is very long and unlikely to have been created by the user before this script is run, it would be replaced
-tempWS = '_ETrans_loading_tempory_WS'
+ monitor1_spec = int(instrument.getNumberParameter("ei-mon1-spec")[0])
+ monitor2_spec = int(instrument.getNumberParameter("ei-mon2-spec")[0])
+
+ # Get incident energy
+ alg = GetEi(input_ws, monitor1_spec, monitor2_spec, E_guess,FixEi=fixei,AdjustBins=True)
+ mon1_peak = float(alg.getPropertyValue("FirstMonitorPeak"))
+ ei = float(input_ws.getSampleDetails().getLogData("Ei").value())
-###########################################
-# Applies unit conversion, detector convcy and grouping correction to a
-# raw file
-###########################################
-def mono_sample(instrum, runs, Ei, d_rebin, wbrf, getEi=True, back='', norma='', det_map='', det_mask='', nameInOut = 'mono_sample_temporyWS') :
- (instru, TOFLow, TOFHigh, run_nums) = retrieveSettings(instrum, back, runs)
+ return ei, mon1_peak
- #----Calculations start------------------
- try:
- pInOut = instru.loadRun(run_nums[0], nameInOut) #this new workspace is accessed by its name nameInOut = 'mono_sample_temporyWS' or pInOut (a pointer)
- if pInOut.isGroup() : raise Exception("Workspace groups are not supported here")
- common.sumWorkspaces(pInOut, instru, run_nums) #the final workspace will take the X-values and instrument from the first workspace and so we don't need to rerun ChangeBinOffset(), MoveInstrumentComponent(), etc.
- Ei, offSet = findPeaksOffSetTimeAndEi(pInOut, Ei, getEi, instru.getFileName(run_nums[0]))
-
- if back != 'noback': #remove the count rate seen in the regions of the histograms defined as the background regions, if the user defined a region
- ConvertToDistribution(pInOut) #deal correctly with changing bin widths
- FlatBackground(pInOut, pInOut, TOFLow-offSet, TOFHigh-offSet, '', 'Mean') #the TOFs in the workspace were offset in a command above so we must take care referring to times now
- ConvertFromDistribution(pInOut) #some algorithms later can't deal with distributions
-
- # deals with normalize to monitor (e.g. norm = 'monitor-monitor 1'), current (if norm = 'protons (uAh)'), etc.
- NormaliseTo(getNorm(instru, norma), pInOut, offSet, instru)
+def NormaliseTo(scheme, data_ws, offset, instrument):
- ConvertUnits(pInOut, pInOut, 'DeltaE', 'Direct', Ei, AlignBins=0)
+ if scheme == 'monitor-monitor 1':
+ min = instrument.getNumberParameter("norm-mon1-min")[0]
+ max = instrument.getNumberParameter("norm-mon1-max")[0]
+ min += offset
+ max += offset
- if d_rebin != '':
- Rebin(pInOut, pInOut, common.listToString(d_rebin))
-
- DetectorEfficiencyCor(pInOut, pInOut, Ei)
-
- if det_mask != '' :
- MaskDetectors(Workspace=pInOut, SpectraList=common.listToString(det_mask))
-
- if det_map != '':
- GroupDetectors(pInOut, pInOut, det_map, KeepUngroupedSpectra=0)
+ mon_spec = int(instrument.getNumberParameter("norm-mon1-spec")[0])
+ NormaliseToMonitor( InputWorkspace=data_ws, OutputWorkspace=data_ws, MonitorSpectrum=mon_spec, IntegrationRangeMin=min, IntegrationRangeMax=max,IncludePartialBins=True)
+ elif scheme == 'monitor-monitor 2':
+ PEAKWIDTH = 4.0/100
+ # the origin of time was set to the location of the peak so to get it's area integrate around the origin, PEAKWIDTH fraction of total time of flight
+ min = (-PEAKWIDTH/2)*20000
+ max = (PEAKWIDTH/2)*20000
+ mon_spec = int(instrument.getNumberParameter("norm-mon2-spec")[0])
+ NormaliseToMonitor( InputWorkspace=data_ws, OutputWorkspace=data_ws, MonitorSpectrum=mon_spec, IntegrationRangeMin=min, IntegrationRangeMax=max,IncludePartialBins=True)
+ elif scheme == 'protons (uAh)':
+ NormaliseByCurrent(InputWorkspace=data_ws, OutputWorkspace=data_ws)
+ elif scheme == 'no normalization':
+ return
+ else:
+ raise RunTimeError('Normalisation scheme ' + reference + ' not found. It must be one of monitor, current, peak or none')
+
+def NormaliseToWhiteBeam(WBRun, mono_ws, mapFile, detMask, scheme, prefix, instrument):
+ wbnorm_name = "_ETrans_norm_tempory_WS"
+ wbnorm_ws = common.LoadNexRaw(common.getFileName(prefix, WBRun), wbnorm_name)[0]
+
+ NormaliseTo(scheme, wbnorm_ws, 0., instrument)
+ ConvertUnits(wbnorm_ws, wbnorm_ws, "Energy", AlignBins=0)
+ # This both integrates the workspace into one bin spectra and sets up common bin boundaries for all spectra
+ low = instrument.getNumberParameter("wb-integr-min")[0]
+ hi = instrument.getNumberParameter("wb-integr-max")[0]
+ delta = 2.0*(hi - low)
+ Rebin(wbnorm_ws, wbnorm_ws, [low, delta, hi])
+
+ if detMask != '' :
+ MaskDetectors(wbnorm_ws, SpectraList=detMask)
+
+ if mapFile!= "" :
+ GroupDetectors(wbnorm_ws, wbnorm_ws, mapFile, KeepUngroupedSpectra=0)
- ConvertToDistribution(pInOut)
- pInOut *= instru.scale_factor
-
- #replaces inifinities and error values with large numbers. Infinity values can be normally be avoided passing good energy values to ConvertUnits
- ReplaceSpecialValues(pInOut, pInOut, 1e40, 1e40, 1e40, 1e40)
+ # White beam scale factor
+ wb_scale_factor = instrument.getNumberParameter("wb-scale-factor")[0]
+ wbnorm_ws *= wb_scale_factor
+ mono_ws /= wbnorm_ws
+
+ mtd.deleteWorkspace(wbnorm_name)
+ return
- NormaliseToWhiteBeam(instru, wbrf, pInOut, det_map, det_mask, norma)
- pInOut /= instru.white_beam_scale
-
- return pInOut
-
- except Exception, reason:
- # delete the possibly part finished workspaces
- for workspace in mantid.getWorkspaceNames() :
- if (workspace == nameInOut) : mantid.deleteWorkspace(nameInOut)
- if (workspace == tempWS) : mantid.deleteWorkspace(tempWS)
- print reason
+# a workspace name that is very long and unlikely to have been created by the user before this script is run, it would be replaced
+tempWS = '_ETrans_loading_tempory_WS'
diff --git a/Code/Mantid/PythonAPI/scripts/Excitations/ExcitDefaults.py b/Code/Mantid/PythonAPI/scripts/Excitations/ExcitDefaults.py
index 6659d534125ffe5383999cdcaaeae9abfd768638..9f4a5953cddd99471caf6b6e0755bd8254391cd1 100644
--- a/Code/Mantid/PythonAPI/scripts/Excitations/ExcitDefaults.py
+++ b/Code/Mantid/PythonAPI/scripts/Excitations/ExcitDefaults.py
@@ -1,26 +1,48 @@
-#! you need to restart the application you are using (e.g. MantidPlot) before changes to this file will take effect!
import CommonFunctions as common
MARI = common.defaults(
- background_range=(18000.0,19500.0),
- normalization='monitor-monitor 1',
- instrument_pref='MAR',
- white_beam_integr=(20.0,40.0),
- scale_factor=1.8182e8,
- monitor1_integr=(1000,2000),
- white_beam_scale=1000.
+ background_range=(18000.0,19500.0),
+ normalization='monitor-monitor 1',
+ instrument_pref='MAR',
+ white_beam_integr=(20.0,40.0),
+ scale_factor=1.8182e8,
+ monitor1_integr=(1000,2000),
+ white_beam_scale=1000.,
+ getei_monitors=(2,3)
)
MAPS = common.defaults(
- background_range=(12000.0, 18000.0),
- normalization='monitor-monitor 1',
- instrument_pref='MAP')
+ background_range=(12000.0,18000.0),
+ normalization='monitor-monitor 1',
+ instrument_pref='MAP',
+ white_beam_integr=(20.0,300.0),
+ scale_factor=1.7016e8,
+ monitor1_integr=(1000,2000),
+ white_beam_scale=1000.,
+ getei_monitors=(41474,41475)
+ )
+
+MERLIN = common.defaults(
+ background_range=(12000.0,18000.0),
+ normalization='monitor-monitor 1',
+ instrument_pref='MER',
+ white_beam_integr=(20.0,55.0),
+ scale_factor=1.7016e8,
+ monitor1_integr=(1000,2000),
+ white_beam_scale=1000.,
+ getei_monitors=(69634,69638)
+ )
-# we only have support for MARI at the moment, ??STEVES?? change this
-def loadDefaults(instPrefix):
- if instPrefix == 'mar' : instPrefix = 'MAR'
- if instPrefix == 'MAR' : return MARI
- else : raise Exception("As only MARI is supported at the moment the first argument must be MAR")
+def loadDefaults(prefix):
+ prefix = prefix.lower()
+ if prefix == 'mar':
+ return MARI
+ elif prefix == 'map':
+ return MAPS
+ elif prefix == 'mer':
+ return MERLIN
+ else:
+ raise RunTimeError("Unknown instrument prefix selected.")
#called r in DIAG
DetectorEfficiencyVariation_Variation = 1.1
\ No newline at end of file
diff --git a/Code/Mantid/PythonAPI/scripts/Excitations/GUI_Interface.py b/Code/Mantid/PythonAPI/scripts/Excitations/GUI_Interface.py
index a87623f5cc59625f7b94eccd174a1d7502119e79..86710e422d3dd7d30f72a0dbc4a3b19d521e5261 100644
--- a/Code/Mantid/PythonAPI/scripts/Excitations/GUI_Interface.py
+++ b/Code/Mantid/PythonAPI/scripts/Excitations/GUI_Interface.py
@@ -1,16 +1,3 @@
-# Used by MantidPlot if the script in is in Code\Mantid\PythonAPI\scripts\Excitations make a
-# copy of the script if you want to change it or risk the MantidPlot interface not working
-###########################################
-# applies unit conversion, detector convcy and grouping correction to a
-# raw file
-###########################################
-#from ConversionLib import *
-#from DetectorTestLib import *
-
-#this commented out code will replace the code immediately below it when the interface been diagnose and the interface is implemented
-#if ??? :
-# runs = '|GUI_SET_RAWFILE_LIST|'
-#badSpectra = diagnose(instrum='MAR',wbrf='11060',wbrf2='11060',runs='15537',tiny=1e-10,huge=1e10,median_lo=0.1,median_hi=3.0,sv_sig=3.3,bmedians=5.0,zero='False', out_asc='',maskFile='')
badSpectra=[]
DetectorMask = ''#+|MASK_WORKSPACE|
if DetectorMask != '' :
@@ -18,5 +5,5 @@ if DetectorMask != '' :
badSpectra = FDOL.getPropertyValue('BadSpectraNums')
mantid.deleteWorkspace('_ETrans_loading_bad_detector_WS')
-mtd.sendLogMessage("--Executing Python function: mono_sample('MAR', '|GUI_SET_RAWFILE_LIST|', |GUI_SET_E_GUESS|, |GUI_SET_BIN_BOUNDS|, |GUI_SET_WBV|, getEi=|GUI_SET_E|, back='('+|TOF_LOW|+','+|TOF_HIGH|+')', norma=|GUI_SET_NORM|, det_map=|GUI_SET_MAP_FILE|, det_mask=badSpectra, nameInOut='|GUI_SET_OUTWS|')--")
-mono_sample('MAR', '|GUI_SET_RAWFILE_LIST|', |GUI_SET_E_GUESS|, |GUI_SET_BIN_BOUNDS|, |GUI_SET_WBV|, getEi=|GUI_SET_E|, back='('+|TOF_LOW|+','+|TOF_HIGH|+')', norma=|GUI_SET_NORM|, det_map=|GUI_SET_MAP_FILE|, det_mask=badSpectra, nameInOut='|GUI_SET_OUTWS|')
+#mtd.sendLogMessage("--Executing Python function: mono_sample('MAR', '|GUI_SET_RAWFILE_LIST|', |GUI_SET_E_GUESS|, |GUI_SET_BIN_BOUNDS|, |GUI_SET_WBV|, getEi=|GUI_SET_E|, back='('+|TOF_LOW|+','+|TOF_HIGH|+')', norma=|GUI_SET_NORM|, det_map=|GUI_SET_MAP_FILE|, det_mask=badSpectra, nameInOut='|GUI_SET_OUTWS|')--")
+mono_sample('|PREFIX|', '|GUI_SET_RAWFILE_LIST|', |GUI_SET_E_GUESS|, |GUI_SET_BIN_BOUNDS|, |GUI_SET_WBV|, getEi=|GUI_SET_E|, back='('+|TOF_LOW|+','+|TOF_HIGH|+')', norma=|GUI_SET_NORM|, det_map=|GUI_SET_MAP_FILE|, det_mask=badSpectra, output_name='|GUI_SET_OUTWS|')
diff --git a/Code/qtiplot/MantidQt/CustomInterfaces/inc/Homer.ui b/Code/qtiplot/MantidQt/CustomInterfaces/inc/Homer.ui
index a9ab9b4eb5804cb9bcc4b2f429812cd49df731b3..94bd06038f3f8596fe76145fe7a6eee61d82e76b 100644
--- a/Code/qtiplot/MantidQt/CustomInterfaces/inc/Homer.ui
+++ b/Code/qtiplot/MantidQt/CustomInterfaces/inc/Homer.ui
@@ -99,7 +99,7 @@
<item row="0" column="0" colspan="2" >
<widget class="QLabel" name="lbPrefix" >
<property name="enabled" >
- <bool>false</bool>
+ <bool>true</bool>
</property>
<property name="text" >
<string>Instrument File Prefix</string>
@@ -109,7 +109,7 @@
<item row="0" column="2" >
<widget class="QComboBox" name="loadRun_cbInst" >
<property name="enabled" >
- <bool>false</bool>
+ <bool>true</bool>
</property>
<property name="sizePolicy" >
<sizepolicy vsizetype="Fixed" hsizetype="Fixed" >
@@ -126,6 +126,21 @@
<property name="editable" >
<bool>true</bool>
</property>
+ <item>
+ <property name="text" >
+ <string>MAR</string>
+ </property>
+ </item>
+ <item>
+ <property name="text" >
+ <string>MAP</string>
+ </property>
+ </item>
+ <item>
+ <property name="text" >
+ <string>MER</string>
+ </property>
+ </item>
</widget>
</item>
<item row="0" column="3" colspan="3" >
diff --git a/Code/qtiplot/MantidQt/CustomInterfaces/src/Homer.cpp b/Code/qtiplot/MantidQt/CustomInterfaces/src/Homer.cpp
index 7595773a4b62ebdd527bf694f5eb0def88ca381b..52885456a8b048b9cf371d134015763159939182 100644
--- a/Code/qtiplot/MantidQt/CustomInterfaces/src/Homer.cpp
+++ b/Code/qtiplot/MantidQt/CustomInterfaces/src/Homer.cpp
@@ -30,10 +30,6 @@ using namespace Mantid::Kernel;
using namespace MantidQt::MantidWidgets;
using namespace MantidQt::CustomInterfaces;
-//these two defaults will be removed when other instruments are supported
-static const QString G_INSTRUMENT("MAR");
-static const QString G_DEFAULT_MAP_FILE("mari_res.map");
-
//default values
static const int G_NUM_NORM_SCHEMES = 3;
static const QString G_NORM_SCHEMES[G_NUM_NORM_SCHEMES] =
@@ -101,22 +97,14 @@ void Homer::pythonIsRunning(bool running)
* and set the current text to the one that was passed
*/
QString Homer::setUpInstru()
-//??STEVES?? move this function to a File widget
-{ // if there were no previously used instruments the "" below adds a blank entry. The empty string entry will always be there, even as more instruments are added
- QStringList prevInstrus =
- m_prev.value("CustomInterfaces/Homer/instrusList","").toStringList();
-
- /*??STEVES ?? get rid of this when more instruments are supported*/if ( ! prevInstrus.contains(G_INSTRUMENT) ) prevInstrus.prepend(G_INSTRUMENT);
- QStringList::const_iterator instru = prevInstrus.begin();
- for ( ; instru != prevInstrus.end(); ++instru )
+{
+ QString curInstru = m_prev.value("CustomInterfaces/Homer/instrument", "").toString();
+ int index = m_uiForm.loadRun_cbInst->findText(curInstru);
+ if( index < 0 )
{
- m_uiForm.loadRun_cbInst->addItem(*instru);
+ index = 0;
}
-
- QString curInstru
- = m_prev.value("CustomInterfaces/Homer/instrument", G_INSTRUMENT).toString();
- m_uiForm.loadRun_cbInst->setEditText(curInstru);
-
+ m_uiForm.loadRun_cbInst->setCurrentIndex(index);
return curInstru;
}
/// For each widgets in the first tab this adds custom widgets, fills in combination boxes and runs setToolTip()
@@ -234,8 +222,7 @@ void Homer::page1Defaults()
m_prev.setValue("TOFend",
m_prev.value("TOFend", G_END_WINDOW_TOF).toDouble());
- m_uiForm.map_fileInput_leName->setText(
- m_prev.value("map", G_DEFAULT_MAP_FILE).toString());
+ //m_uiForm.map_fileInput_leName->setText(m_prev.value("map", G_DEFAULT_MAP_FILE).toString());
}
/// make validator labels and associate them with the controls that need them in the first tab
void Homer::page1Validators()
diff --git a/Code/qtiplot/MantidQt/CustomInterfaces/src/deltaECalc.cpp b/Code/qtiplot/MantidQt/CustomInterfaces/src/deltaECalc.cpp
index 913a685d023e5d38a3aa5fc91827500b602a7309..bee077848b1e7f9694e03828db123baadbd83314 100644
--- a/Code/qtiplot/MantidQt/CustomInterfaces/src/deltaECalc.cpp
+++ b/Code/qtiplot/MantidQt/CustomInterfaces/src/deltaECalc.cpp
@@ -102,6 +102,7 @@ QString deltaECalc::createProcessingScript(const std::string &inFiles, const std
QString err;
+ m_pyScript.replace("|PREFIX|", m_sets.loadRun_cbInst->currentText());
// here we are placing code directly into specified parts of the Python that will get run
m_pyScript.replace("|GUI_SET_RAWFILE_LIST|", QString::fromStdString(inFiles));
// these functions replace whole blocks of code
@@ -126,6 +127,8 @@ QString deltaECalc::createProcessingScript(const std::string &inFiles, const std
m_pyScript.replace("|TEMPWS|", tempWS);
+ std::cerr << m_pyScript.toStdString() << "\n";
+
return WSName;
}
/** Completes the Python statements that implement calculating the initial energy