Skip to content
Snippets Groups Projects
Commit afbcfd9b authored by Zhou, Wenduo's avatar Zhou, Wenduo
Browse files

Implementing FilterEvents and etc. Refs #5056.

1. Add algorithm FilterEvents to filter events from SplittersWorkspace
2. Add but not finish algorithm GenerateEventsFilter to support
FilterEvents
3. Enhance EventList, SplittersWorkspace, WorkspaceCreationHelper to
support the new algorithms.
parent bd0c7d02
No related merge requests found
Showing
with 1272 additions and 12 deletions
...@@ -79,6 +79,7 @@ set ( SRC_FILES ...@@ -79,6 +79,7 @@ set ( SRC_FILES
src/FilterBadPulses.cpp src/FilterBadPulses.cpp
src/FilterByLogValue.cpp src/FilterByLogValue.cpp
src/FilterByTime.cpp src/FilterByTime.cpp
src/FilterEvents.cpp
src/FilterEventsHighFrequency.cpp src/FilterEventsHighFrequency.cpp
src/FindCenterOfMassPosition.cpp src/FindCenterOfMassPosition.cpp
src/FindCenterOfMassPosition2.cpp src/FindCenterOfMassPosition2.cpp
...@@ -88,6 +89,8 @@ set ( SRC_FILES ...@@ -88,6 +89,8 @@ set ( SRC_FILES
src/FlatBackground.cpp src/FlatBackground.cpp
src/FlatPlateAbsorption.cpp src/FlatPlateAbsorption.cpp
src/GeneralisedSecondDifference.cpp src/GeneralisedSecondDifference.cpp
src/GenerateEventsFilter.cpp
src/GeneratePeaks.cpp
src/GeneratePythonScript.cpp src/GeneratePythonScript.cpp
src/GetDetOffsetsMultiPeaks.cpp src/GetDetOffsetsMultiPeaks.cpp
src/GetDetectorOffsets.cpp src/GetDetectorOffsets.cpp
...@@ -267,6 +270,7 @@ set ( INC_FILES ...@@ -267,6 +270,7 @@ set ( INC_FILES
inc/MantidAlgorithms/FilterBadPulses.h inc/MantidAlgorithms/FilterBadPulses.h
inc/MantidAlgorithms/FilterByLogValue.h inc/MantidAlgorithms/FilterByLogValue.h
inc/MantidAlgorithms/FilterByTime.h inc/MantidAlgorithms/FilterByTime.h
inc/MantidAlgorithms/FilterEvents.h
inc/MantidAlgorithms/FilterEventsHighFrequency.h inc/MantidAlgorithms/FilterEventsHighFrequency.h
inc/MantidAlgorithms/FindCenterOfMassPosition.h inc/MantidAlgorithms/FindCenterOfMassPosition.h
inc/MantidAlgorithms/FindCenterOfMassPosition2.h inc/MantidAlgorithms/FindCenterOfMassPosition2.h
...@@ -277,6 +281,8 @@ set ( INC_FILES ...@@ -277,6 +281,8 @@ set ( INC_FILES
inc/MantidAlgorithms/FlatPlateAbsorption.h inc/MantidAlgorithms/FlatPlateAbsorption.h
inc/MantidAlgorithms/GSLFunctions.h inc/MantidAlgorithms/GSLFunctions.h
inc/MantidAlgorithms/GeneralisedSecondDifference.h inc/MantidAlgorithms/GeneralisedSecondDifference.h
inc/MantidAlgorithms/GenerateEventsFilter.h
inc/MantidAlgorithms/GeneratePeaks.h
inc/MantidAlgorithms/GeneratePythonScript.h inc/MantidAlgorithms/GeneratePythonScript.h
inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h inc/MantidAlgorithms/GetDetOffsetsMultiPeaks.h
inc/MantidAlgorithms/GetDetectorOffsets.h inc/MantidAlgorithms/GetDetectorOffsets.h
...@@ -376,6 +382,7 @@ set ( INC_FILES ...@@ -376,6 +382,7 @@ set ( INC_FILES
) )
set ( TEST_FILES set ( TEST_FILES
# test/CalMuonDeadTimeTest.h
test/AddLogDerivativeTest.h test/AddLogDerivativeTest.h
test/AddSampleLogTest.h test/AddSampleLogTest.h
test/AlignDetectorInTOFTest.h test/AlignDetectorInTOFTest.h
...@@ -389,7 +396,6 @@ set ( TEST_FILES ...@@ -389,7 +396,6 @@ set ( TEST_FILES
test/BinaryOperateMasksTest.h test/BinaryOperateMasksTest.h
test/BinaryOperationTest.h test/BinaryOperationTest.h
test/BlendSqTest.h test/BlendSqTest.h
# test/CalMuonDeadTimeTest.h
test/CalculateEfficiencyTest.h test/CalculateEfficiencyTest.h
test/CalculateTransmissionBeamSpreaderTest.h test/CalculateTransmissionBeamSpreaderTest.h
test/CalculateTransmissionTest.h test/CalculateTransmissionTest.h
...@@ -449,6 +455,7 @@ set ( TEST_FILES ...@@ -449,6 +455,7 @@ set ( TEST_FILES
test/FilterByLogValueTest.h test/FilterByLogValueTest.h
test/FilterByTimeTest.h test/FilterByTimeTest.h
test/FilterEventsHighFrequencyTest.h test/FilterEventsHighFrequencyTest.h
test/FilterEventsTest.h
test/FindCenterOfMassPosition2Test.h test/FindCenterOfMassPosition2Test.h
test/FindCenterOfMassPositionTest.h test/FindCenterOfMassPositionTest.h
test/FindDeadDetectorsTest.h test/FindDeadDetectorsTest.h
...@@ -456,6 +463,8 @@ set ( TEST_FILES ...@@ -456,6 +463,8 @@ set ( TEST_FILES
test/FindPeaksTest.h test/FindPeaksTest.h
test/FlatBackgroundTest.h test/FlatBackgroundTest.h
test/FlatPlateAbsorptionTest.h test/FlatPlateAbsorptionTest.h
test/GenerateEventsFilterTest.h
test/GeneratePeaksTest.h
test/GeneratePythonScriptTest.h test/GeneratePythonScriptTest.h
test/GetDetOffsetsMultiPeaksTest.h test/GetDetOffsetsMultiPeaksTest.h
test/GetDetectorOffsetsTest.h test/GetDetectorOffsetsTest.h
......
#ifndef MANTID_ALGORITHMS_FILTEREVENTS_H_
#define MANTID_ALGORITHMS_FILTEREVENTS_H_
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/SplittersWorkspace.h"
#include "MantidAPI/ISplittersWorkspace.h"
#include "MantidKernel/TimeSplitter.h"
namespace Mantid
{
namespace Algorithms
{
/** FilterEvents : Filter Events in EventWorkspace to multiple EventsWorkspace by Splitters
@date 2012-04-04
Copyright © 2012 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
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://svn.mantidproject.org/mantid/trunk/Code/Mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport FilterEvents : public API::Algorithm
{
public:
FilterEvents();
virtual ~FilterEvents();
/// Algorithm's name for identification overriding a virtual method
virtual const std::string name() const { return "FilterEvents";};
/// Algorithm's version for identification overriding a virtual method
virtual int version() const { return 1;};
/// Algorithm's category for identification overriding a virtual method
virtual const std::string category() const { return "Events\\EventFiltering";}
private:
/// Sets documentation strings for this algorithm
virtual void initDocs();
// Implement abstract Algorithm methods
void init();
// Implement abstract Algorithm methods
void exec();
void processSplittersWorkspace();
void createOutputWorkspaces(std::string outputwsnamebase);
void importDetectorTOFCalibration(std::string detcalfilename);
void filterEventsBySplitters();
DataObjects::EventWorkspace_sptr mEventWorkspace;
DataObjects::SplittersWorkspace_sptr mSplittersWorkspace;
std::set<int> mWorkspaceGroups;
Kernel::TimeSplitterType mSplitters;
std::map<int, DataObjects::EventWorkspace_sptr> mOutputWorkspaces;
std::vector<detid_t> mCalibDetectorIDs;
std::vector<double> mCalibOffsets;
bool mFilterByPulseTime;
};
} // namespace Algorithms
} // namespace Mantid
#endif /* MANTID_ALGORITHMS_FILTEREVENTS_H_ */
#ifndef MANTID_ALGORITHMS_GENERATEEVENTSFILTER_H_
#define MANTID_ALGORITHMS_GENERATEEVENTSFILTER_H_
#include "MantidKernel/System.h"
#include "MantidAPI/Algorithm.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidDataObjects/SplittersWorkspace.h"
namespace Mantid
{
namespace Algorithms
{
/** GenerateEventsFilter : Generate an events-filter, i.e., a SplittersWorkspace according
to user's request.
Request can be a combination of log value and time, i.e.,
(1) T_min
(2) T_max
(3) delta Time
(4) Min log value
(5) Max log value
(6) delta log value
(7) identify log value increment (bool)
(8) number of sections per interval (applied to log value only!)
This algorithm can generate filters including
(1) deltaT per interval from T_min to T_max with : Log value is not given
(2) delta log value per interval from T_min to T_max and from min log value to max log value
Note:
(1) Time can be (a) relative time in ns (b) relative time in second (float) (c) percentage time
(2) if option "identify log value increment"
@date 2012-04-09
Copyright &copy; 2012 ISIS Rutherford Appleton Laboratory & NScD Oak Ridge National Laboratory
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://svn.mantidproject.org/mantid/trunk/Code/Mantid>
Code Documentation is available at: <http://doxygen.mantidproject.org>
*/
class DLLExport GenerateEventsFilter : public API::Algorithm
{
public:
GenerateEventsFilter();
virtual ~GenerateEventsFilter();
/// Algorithm's name for identification overriding a virtual method
virtual const std::string name() const { return "GenerateEventsFilter";};
/// Algorithm's version for identification overriding a virtual method
virtual int version() const { return 1;};
/// Algorithm's category for identification overriding a virtual method
virtual const std::string category() const { return "Events\\EventFiltering";}
private:
/// Sets documentation strings for this algorithm
virtual void initDocs();
// Implement abstract Algorithm methods
void init();
// Implement abstract Algorithm methods
void exec();
void processInputTime(Kernel::DateAndTime runstarttime);
void setFilterByTimeOnly();
void setFilterByValue();
DataObjects::EventWorkspace_const_sptr mEventWS;
DataObjects::SplittersWorkspace_sptr mSplitters;
Kernel::DateAndTime mStartTime;
Kernel::DateAndTime mStopTime;
};
} // namespace Algorithms
} // namespace Mantid
#endif /* MANTID_ALGORITHMS_GENERATEEVENTSFILTER_H_ */
#include "MantidAlgorithms/FilterEvents.h"
#include "MantidKernel/System.h"
#include "MantidAPI/WorkspaceProperty.h"
#include "MantidAPI/FileProperty.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidAPI/WorkspaceFactory.h"
#include <sstream>
using namespace Mantid::Kernel;
using namespace Mantid::API;
namespace Mantid
{
namespace Algorithms
{
DECLARE_ALGORITHM(FilterEvents)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
FilterEvents::FilterEvents()
{
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
FilterEvents::~FilterEvents()
{
}
void FilterEvents::initDocs()
{
}
/*
* Declare Inputs
*/
void FilterEvents::init()
{
declareProperty(
new API::WorkspaceProperty<DataObjects::EventWorkspace>("InputWorkspace","",Direction::Input),
"An input event workspace" );
declareProperty("OutputWorkspaceBaseName", "OutputWorkspace",
"The base name to use for the output workspace" );
declareProperty(
new API::WorkspaceProperty<DataObjects::SplittersWorkspace>("InputSplittersWorkspace", "", Direction::Input),
"An input SpilltersWorskpace for filtering");
this->declareProperty(new API::FileProperty("DetectorCalibrationFile", "", API::FileProperty::OptionalLoad, ".dat"),
"Input pixel TOF calibration file in column data format");
this->declareProperty("FilterByPulseTime", false,
"Filter the event by its pulse time only. This option can make execution of algorithm faster. But it lowers precision.");
return;
}
/*
* Execution body
*/
void FilterEvents::exec()
{
// 1. Get inputs
mEventWorkspace = this->getProperty("InputWorkspace");
mSplittersWorkspace = this->getProperty("InputSplittersWorkspace");
std::string outputwsnamebase = this->getProperty("OutputWorkspaceBaseName");
std::string detcalfilename = this->getProperty("DetectorCalibrationFile");
mFilterByPulseTime = this->getProperty("FilterByPulseTime");
// 2. Process inputs
processSplittersWorkspace();
createOutputWorkspaces(outputwsnamebase);
importDetectorTOFCalibration(detcalfilename);
// 3. Filter Events
filterEventsBySplitters();
return;
}
/*
* Convert SplitterWorkspace object to TimeSplitterType (sorted vector)
* and create a map for all workspace group number
*/
void FilterEvents::processSplittersWorkspace()
{
// 1. Init data structure
size_t numsplitters = mSplittersWorkspace->getNumberSplitters();
mSplitters.reserve(numsplitters);
// 2. Insert all splitters
bool inorder = true;
for (size_t i = 0; i < numsplitters; i ++)
{
mSplitters.push_back(mSplittersWorkspace->getSplitter(i));
mWorkspaceGroups.insert(mSplitters.back().index());
if (inorder && i > 0 && mSplitters[i] < mSplitters[i-1])
inorder = false;
}
// 3. Order if not ordered
if (!inorder)
{
std::sort(mSplitters.begin(), mSplitters.end());
}
mWorkspaceGroups.insert(-1);
return;
}
/*
* Create a list of EventWorkspace for output
*/
void FilterEvents::createOutputWorkspaces(std::string outputwsnamebase)
{
std::set<int>::iterator groupit;
for (groupit = mWorkspaceGroups.begin(); groupit != mWorkspaceGroups.end(); ++groupit)
{
// 1. Get workspace name
int wsgroup = *groupit;
std::stringstream wsname;
wsname << outputwsnamebase << "_" << wsgroup;
std::stringstream parname;
parname << "OutputWorkspace_" << wsgroup;
// 2. Generate one of the output workspaces & Copy geometry over. But we don't copy the data.
DataObjects::EventWorkspace_sptr optws = boost::dynamic_pointer_cast<DataObjects::EventWorkspace>(
API::WorkspaceFactory::Instance().create("EventWorkspace", mEventWorkspace->getNumberHistograms(), 2, 1));
API::WorkspaceFactory::Instance().initializeFromParent(mEventWorkspace, optws, false);
// 3. Set to map
mOutputWorkspaces.insert(std::make_pair(wsgroup, optws));
// 4. Set to output workspace
this->declareProperty(new API::WorkspaceProperty<DataObjects::EventWorkspace>(parname.str(), wsname.str(), Direction::Output), "Output");
this->setProperty(parname.str(), optws);
std::cout << "DB9141 Output Workspace: Group = " << wsgroup << " Property Name = " << wsname.str() << std::endl;
} // ENDFOR
return;
}
/*
* Import the detector calibration on TOF
*/
void FilterEvents::importDetectorTOFCalibration(std::string detcalfilename)
{
detid_t indet;
// 1. Check workspace
if (!mEventWorkspace)
{
g_log.error() << "Required to import EventWorkspace before calling importCalibrationFile()" << std::endl;
throw std::invalid_argument("Calling function in wrong order!");
}
// 2. Prepare output
mCalibDetectorIDs.clear();
mCalibOffsets.clear();
size_t numhist = mEventWorkspace->getNumberHistograms();
mCalibDetectorIDs.reserve(numhist);
mCalibOffsets.reserve(numhist);
// 3. Read file?
bool readcalfile = true;
if (detcalfilename.empty())
{
readcalfile = false;
}
if (readcalfile)
{
try{
// a. Open file
std::ifstream ifs;
ifs.open(detcalfilename.c_str(), std::ios::in);
double doffset;
for (size_t i = 0; i < numhist; i ++)
{
// i. each pixel: get detector ID from EventWorkspace
const DataObjects::EventList events = mEventWorkspace->getEventList(i);
std::set<detid_t> detids = events.getDetectorIDs();
std::set<detid_t>::iterator detit;
detid_t detid = 0;
for (detit=detids.begin(); detit!=detids.end(); ++detit)
detid = *detit;
// ii. read file
ifs >> indet >> doffset;
// iii. store
if (indet != detid){
g_log.error() << "Calibration File Error! Line " << i << " should read in pixel " << detid << " but read in " << indet
<< "\nAbort to reading calibration file!"<< std::endl;
readcalfile = false;
break;
}
else if (doffset < 0 || doffset > 1.0)
{
g_log.error() << "Calibration File Error! Line " << i << " have ratio offset outside (0,1) " << detid << " but read in " << indet
<<"\nAbort to reading calibration file!"<< std::endl;
readcalfile = false;
break;
}
else
{
mCalibDetectorIDs.push_back(detid);
mCalibOffsets.push_back(doffset);
}
}
ifs.close();
}
catch (std::ifstream::failure&)
{
g_log.error() << "Calibration File Error! Open calibration/offset file " << detcalfilename << " error " << std::endl;
mCalibDetectorIDs.clear();
mCalibOffsets.clear();
readcalfile = false;
}
} // If-readcalfile
// 4. Use default/dummy offset calibration = 1.0
if (!readcalfile)
{
g_log.notice() << "Using default detector offset/calibration" << std::endl;
for (size_t i = 0; i < mEventWorkspace->getNumberHistograms(); i ++)
{
std::set<detid_t> detids = mEventWorkspace->getEventList(i).getDetectorIDs();
std::set<detid_t>::iterator detit;
detid_t detid = 0;
for (detit=detids.begin(); detit!=detids.end(); ++detit)
detid = *detit;
mCalibDetectorIDs.push_back(detid);
mCalibOffsets.push_back(1.0);
}
} // If NOT Read-calibration-file
return;
}
/*
* Main filtering method
*/
void FilterEvents::filterEventsBySplitters()
{
size_t numberOfSpectra = mEventWorkspace->getNumberHistograms();
std::map<int, DataObjects::EventWorkspace_sptr>::iterator wsiter;
// 1. Loop over the histograms (detector spectra)
// PARALLEL_FOR_NO_WSP_CHECK()
for (int64_t i = 0; i < int64_t(numberOfSpectra); ++i)
{
// PARALLEL_START_INTERUPT_REGION
// a) Get the output event lists (should be empty) to be a map
std::map<int, DataObjects::EventList* > outputs;
for (wsiter = mOutputWorkspaces.begin(); wsiter != mOutputWorkspaces.end(); ++ wsiter)
{
int index = wsiter->first;
DataObjects::EventList* output_el = wsiter->second->getEventListPtr(i);
outputs.insert(std::make_pair(index, output_el));
}
// b) and this is the input event list
const DataObjects::EventList& input_el = mEventWorkspace->getEventList(i);
// c) Perform the filtering (using the splitting function and just one output)
input_el.splitByFullTime(mSplitters, outputs, mCalibOffsets[i]);
// prog.report();
// PARALLEL_END_INTERUPT_REGION
}
// PARALLEL_CHECK_INTERUPT_REGION
// 2. Finish adding events and To split/filter the runs,
for (wsiter = mOutputWorkspaces.begin(); wsiter != mOutputWorkspaces.end(); ++wsiter)
{
// 2a Done adding event
wsiter->second->doneAddingEventLists();
// 2b To split/filter the runs, make a vector with just the one output run
// FIXME Complete split the "Run".
std::vector< Run *> output_runs;
output_runs.push_back( &wsiter->second->mutableRun() );
mEventWorkspace->run().splitByTime(mSplitters, output_runs);
}
return;
}
} // namespace Mantid
} // namespace Algorithms
#include "MantidAlgorithms/GenerateEventsFilter.h"
#include "MantidKernel/System.h"
#include "MantidKernel/ListValidator.h"
using namespace Mantid::Kernel;
using namespace Mantid::API;
namespace Mantid
{
namespace Algorithms
{
DECLARE_ALGORITHM(GenerateEventsFilter)
//----------------------------------------------------------------------------------------------
/** Constructor
*/
GenerateEventsFilter::GenerateEventsFilter()
{
}
//----------------------------------------------------------------------------------------------
/** Destructor
*/
GenerateEventsFilter::~GenerateEventsFilter()
{
}
void GenerateEventsFilter::initDocs()
{
}
/*
* Define input
*/
void GenerateEventsFilter::init()
{
declareProperty(
new WorkspaceProperty<DataObjects::EventWorkspace>("InputWorkspace","",Direction::InOut),
"An input event workspace" );
declareProperty(
new WorkspaceProperty<DataObjects::SplittersWorkspace>("OutputWorkspace","",Direction::Output),
"The name to use for the output SplittersWorkspace object, i.e., the filter." );
// 1. Time
declareProperty("StartTime", 0.0,
"The start time, in (a) seconds, (b) nanoseconds or (c) percentage of total run time\n"
"since the start of the run. Events before this time are filtered out. \n"
"run_start is used as the zero. ");
declareProperty("StopTime", -1.0,
"The stop time, in (2) seconds, (b) nanoseconds or (c) percentage of total run time\n"
"since the start of the run. Events at or after this time are filtered out. \n"
"run_start is used as the zero. ");
declareProperty("TimeInterval", -1.0,
"Length of the time splices if filtered in time only.");
std::vector<std::string> timeoptions;
timeoptions.push_back("Seconds");
timeoptions.push_back("Nanoseconds");
timeoptions.push_back("Percent");
declareProperty("TimeType", "Seconds", boost::make_shared<Kernel::StringListValidator>(timeoptions),
"Type to define StartTime, StopTime and DeltaTime");
// 2. Log value
declareProperty("LogName", "",
"Name of the sample log to use to filter.\n"
"For example, the pulse charge is recorded in 'ProtonCharge'.");
declareProperty("MinimumValue", 0.0, "Minimum log value for which to keep events.");
declareProperty("MaximumValue", -1.0, "Maximum log value for which to keep events.");
declareProperty("DeltaValue", -1.0,
"Delta of log value to be sliced into from min log value and max log value.\n"
"If not given, then only value ");
declareProperty("SeparateLogValueChangeDirection", false,
"d(log value)/dt can be positive and negative. They can be put to different splitters.");
declareProperty("LogValueTimeSections", 1,
"In one log value interval, it can be further divided into sections in even time slice.");
}
/*
* Main execute body
*/
void GenerateEventsFilter::exec()
{
// 1. Get general input and output
mEventWS = this->getProperty("InputWorkspace");
Kernel::DateAndTime runstart(mEventWS->run().getProperty("run_start")->value());
mSplitters = boost::shared_ptr<DataObjects::SplittersWorkspace>(new DataObjects::SplittersWorkspace);
this->setProperty("OutputWorkspace", mSplitters);
std::cout << "DB9441 Run Start = " << runstart << " / " << runstart.totalNanoseconds() << std::endl;
// 2. Get Time
processInputTime(runstart);
// 3. Get Log
std::string logname = this->getProperty("LogName");
if (logname.empty())
{
// a) Set filter by time only
setFilterByTimeOnly();
}
else
{
// b) Set filter by time and log
setFilterByValue();
}
}
/*
* Process the input for time. A smart but complicated default rule
*/
void GenerateEventsFilter::processInputTime(Kernel::DateAndTime runstarttime)
{
// 1. Get input
double inpt0 = this->getProperty("StartTime");
double inptf = this->getProperty("StopTime");
std::string timetype = this->getProperty("TimeType");
if (inpt0 < 0)
{
throw std::invalid_argument("Input StartTime cannot be negative!");
}
// 2. Find maximum time by proton charge
// FIXME Use this simple method may miss the events in the last pulse
Kernel::TimeSeriesProperty<double>* protonchargelog =
dynamic_cast<Kernel::TimeSeriesProperty<double> *>(mEventWS->run().getProperty("proton_charge"));
Kernel::DateAndTime runend = protonchargelog->lastTime();
// 3. Set up time-convert unit
double convertfactor = 1.0;
if (timetype.compare("Seconds") == 0)
{
// a) In unit of seconds
convertfactor = 1.0E9;
}
else if (timetype.compare("Nanoseconds") == 0)
{
// b) In unit of nano-seconds
convertfactor = 1.0;
}
else if (timetype.compare("Percent") == 0)
{
// c) In unit of percent of total run time
int64_t runtime_ns = runend.totalNanoseconds()-runstarttime.totalNanoseconds();
double runtimed_ns = static_cast<double>(runtime_ns);
convertfactor = 0.01*runtimed_ns;
}
else
{
// d) Not defined
g_log.error() << "TimeType " << timetype << " is not supported!" << std::endl;
throw std::runtime_error("Input TimeType is not supported");
}
// 4. Process second round
int64_t t0_ns = runstarttime.totalNanoseconds() + static_cast<int64_t>(inpt0*convertfactor);
Kernel::DateAndTime tmpt0(t0_ns);
mStartTime = tmpt0;
if (inptf < inpt0)
{
mStopTime = runend;
}
else
{
int64_t tf_ns = runstarttime.totalNanoseconds()+static_cast<int64_t>(inptf*convertfactor);
Kernel::DateAndTime tmptf(tf_ns);
mStopTime = tmptf;
}
return;
}
/*
* Set splitters by time value / interval only
*/
void GenerateEventsFilter::setFilterByTimeOnly()
{
double timeinterval = this->getProperty("TimeInterval");
if (timeinterval <= 0.0)
{
// 1. Default and thus just one interval
}
}
void GenerateEventsFilter::setFilterByValue()
{
// Kernel::TimeSeriesProperty<double>* mLog =
}
} // namespace Mantid
} // namespace Algorithms
#ifndef MANTID_ALGORITHMS_FILTEREVENTSTEST_H_
#define MANTID_ALGORITHMS_FILTEREVENTSTEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
#include <iostream>
#include <iomanip>
#include "MantidAlgorithms/FilterEvents.h"
#include "MantidTestHelpers/WorkspaceCreationHelper.h"
#include "MantidDataObjects/EventWorkspace.h"
#include "MantidDataObjects/Events.h"
#include "MantidDataObjects/EventList.h"
#include "MantidDataObjects/SplittersWorkspace.h"
#include "MantidKernel/TimeSplitter.h"
using namespace Mantid;
using namespace Mantid::Algorithms;
using namespace Mantid::API;
class FilterEventsTest : 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 FilterEventsTest *createSuite() { return new FilterEventsTest(); }
static void destroySuite( FilterEventsTest *suite ) { delete suite; }
void test_Initialization()
{
FilterEvents alg;
alg.initialize();
TS_ASSERT(alg.isInitialized());
}
void test_CreatedEventWorskpaceAndSplitter()
{
DataObjects::EventWorkspace_sptr eventws = createEventWorkspace();
eventws->setName("Test01");
TS_ASSERT_EQUALS(eventws->getNumberEvents(), 500);
DataObjects::EventList elist = eventws->getEventList(0);
TS_ASSERT_EQUALS(elist.getNumberEvents(), 50);
TS_ASSERT(elist.hasDetectorID(1));
DataObjects::SplittersWorkspace_sptr splittersws = createSplitter();
TS_ASSERT_EQUALS(splittersws->getNumberSplitters(), 5);
return;
}
/*
* Filter events without any correction
* (1) Leave correction file empty
* (2) Count events in each output including "-1", the excluded/unselected events
*/
void test_FilterWOCorrection()
{
// 1. Create EventWorkspace and SplittersWorkspace
DataObjects::EventWorkspace_sptr inpWS = createEventWorkspace();
inpWS->setName("Test02");
inpWS->setTitle("Test02");
DataObjects::SplittersWorkspace_sptr splws = createSplitter();
splws->setName("Splitter02");
splws->setTitle("Splitter02");
FilterEvents filter;
filter.initialize();
// 2. Set properties
filter.setProperty("InputWorkspace", inpWS);
filter.setProperty("OutputWorkspaceBaseName", "FilteredWS01");
filter.setProperty("InputSplittersWorkspace", splws);
filter.setProperty("DetectorCalibrationFile", "");
// 3. Execute
TS_ASSERT_THROWS_NOTHING(filter.execute());
TS_ASSERT(filter.isExecuted());
// 4. Get output
// 4.1 Workspace group 0
DataObjects::EventWorkspace_sptr filteredws0 = boost::dynamic_pointer_cast
<DataObjects::EventWorkspace>(AnalysisDataService::Instance().retrieve("FilteredWS01_0"));
TS_ASSERT(filteredws0);
TS_ASSERT_EQUALS(filteredws0->getNumberHistograms(), 10);
TS_ASSERT_EQUALS(filteredws0->getEventList(0).getNumberEvents(), 4);
// 4.2 Workspace group 1
DataObjects::EventWorkspace_sptr filteredws1 = boost::dynamic_pointer_cast
<DataObjects::EventWorkspace>(AnalysisDataService::Instance().retrieve("FilteredWS01_1"));
TS_ASSERT(filteredws1);
TS_ASSERT_EQUALS(filteredws1->getEventList(1).getNumberEvents(), 16);
// 4.3 Workspace group 2
DataObjects::EventWorkspace_sptr filteredws2 = boost::dynamic_pointer_cast
<DataObjects::EventWorkspace>(AnalysisDataService::Instance().retrieve("FilteredWS01_2"));
TS_ASSERT(filteredws2);
TS_ASSERT_EQUALS(filteredws2->getEventList(1).getNumberEvents(), 21);
int64_t runstart_i64 = 20000000000;
int64_t pulsedt = 100*1000*1000;
int64_t tofdt = 10*1000*1000;
DataObjects::EventList elist3 = filteredws2->getEventList(3);
elist3.sortPulseTimeTOF();
DataObjects::TofEvent eventmin = elist3.getEvent(0);
TS_ASSERT_EQUALS(eventmin.pulseTime().totalNanoseconds(), runstart_i64+pulsedt*2);
TS_ASSERT_DELTA(eventmin.tof(), 0, 1.0E-4);
DataObjects::TofEvent eventmax = elist3.getEvent(20);
TS_ASSERT_EQUALS(eventmax.pulseTime().totalNanoseconds(), runstart_i64+pulsedt*4);
TS_ASSERT_DELTA(eventmax.tof(), static_cast<double>(tofdt*6/1000), 1.0E-4);
return;
}
void test_FilterWithCorrection()
{
}
/*
* Create an EventWorkspace. This workspace will have
* (1) Events with wall time even in time
*/
DataObjects::EventWorkspace_sptr createEventWorkspace()
{
// 1. Create an EventWorkspace with 10 detectors
DataObjects::EventWorkspace_sptr eventWS = WorkspaceCreationHelper::createEventWorkspaceWithFullInstrument(10, 1, true);
int64_t runstart_i64 = 20000000000;
Kernel::DateAndTime runstart(runstart_i64);
int64_t pulsedt = 100*1000*1000;
int64_t tofdt = 10*1000*1000;
// 2. Set run_start time
eventWS->mutableRun().addProperty("run_start", runstart.toISO8601String(), true);
for (size_t i = 0; i < eventWS->getNumberHistograms(); i ++)
{
DataObjects::EventList* elist = eventWS->getEventListPtr(i);
for (int64_t pid = 0; pid < 5; pid ++)
{
int64_t pulsetime_i64 = pid*pulsedt+runstart.totalNanoseconds();
Kernel::DateAndTime pulsetime(pulsetime_i64);
for (size_t e = 0; e < 10; e ++)
{
double tof = static_cast<double>(e*tofdt/1000);
DataObjects::TofEvent event(tof, pulsetime);
elist->addEventQuickly(event);
}
} // FOR each pulse
} // For each bank
return eventWS;
}
/*
* Create a Splitter for output
* Region:
* 0: pulse 0: 0 ~ 3+
* 1: pulse 0: 3+ ~ pulse 1: 9+
* 2: from pulse 2: 0 ~ 6+
* -1: from pulse 2: 6+ ~ 9+
*/
DataObjects::SplittersWorkspace_sptr createSplitter()
{
DataObjects::SplittersWorkspace_sptr splitterws = boost::shared_ptr<DataObjects::SplittersWorkspace>(new DataObjects::SplittersWorkspace);
int64_t runstart_i64 = 20000000000;
int64_t pulsedt = 100*1000*1000;
int64_t tofdt = 10*1000*1000;
// 1. Splitter 0: 0 ~ 3+ (first pulse)
int64_t t0 = runstart_i64;
int64_t t1 = t0 + tofdt*3 + tofdt/2;
Kernel::SplittingInterval interval0(t0, t1, 0);
splitterws->addSplitter(interval0);
// 2. Splitter 1: 3+ ~ 9+ (second pulse)
t0 = t1;
t1 = runstart_i64 + pulsedt + tofdt*9 + tofdt/2;
Kernel::SplittingInterval interval1(t0, t1, 1);
splitterws->addSplitter(interval1);
// 3. Splitter 2: from 3rd pulse, 0 ~ 6+
for (size_t i = 2; i < 5; i ++)
{
t0 = runstart_i64+i*pulsedt;
t1 = runstart_i64+i*pulsedt+6*tofdt+tofdt/2;
Kernel::SplittingInterval interval2(t0, t1, 2);
splitterws->addSplitter(interval2);
}
return splitterws;
}
};
#endif /* MANTID_ALGORITHMS_FILTEREVENTSTEST_H_ */
#ifndef MANTID_ALGORITHMS_GENERATEEVENTSFILTERTEST_H_
#define MANTID_ALGORITHMS_GENERATEEVENTSFILTERTEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidKernel/Timer.h"
#include "MantidKernel/System.h"
#include <iostream>
#include <iomanip>
#include "MantidAlgorithms/GenerateEventsFilter.h"
using namespace Mantid;
using namespace Mantid::Algorithms;
using namespace Mantid::API;
class GenerateEventsFilterTest : 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 GenerateEventsFilterTest *createSuite() { return new GenerateEventsFilterTest(); }
static void destroySuite( GenerateEventsFilterTest *suite ) { delete suite; }
void test_Init()
{
GenerateEventsFilter alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize());
TS_ASSERT(alg.isInitialized());
}
};
#endif /* MANTID_ALGORITHMS_GENERATEEVENTSFILTERTEST_H_ */
...@@ -357,6 +357,11 @@ public: ...@@ -357,6 +357,11 @@ public:
void splitByTimeHelper(Kernel::TimeSplitterType & splitter, std::vector< EventList * > outputs, typename std::vector<T> & events) const; void splitByTimeHelper(Kernel::TimeSplitterType & splitter, std::vector< EventList * > outputs, typename std::vector<T> & events) const;
void splitByTime(Kernel::TimeSplitterType & splitter, std::vector< EventList * > outputs) const; void splitByTime(Kernel::TimeSplitterType & splitter, std::vector< EventList * > outputs) const;
template< class T >
void splitByFullTimeHelper(Kernel::TimeSplitterType & splitter, std::map<int, EventList * > outputs, typename std::vector<T> & events,
double tofcorrection) const;
void splitByFullTime(Kernel::TimeSplitterType & splitter, std::map<int, EventList * > outputs, double tofcorrection) const;
template< class T> template< class T>
static void multiplyHelper(std::vector<T> & events, const double value, const double error = 0.0); static void multiplyHelper(std::vector<T> & events, const double value, const double error = 0.0);
void multiply(const double value, const double error = 0.0); void multiply(const double value, const double error = 0.0);
...@@ -381,6 +386,9 @@ public: ...@@ -381,6 +386,9 @@ public:
void convertUnitsViaTof(Mantid::Kernel::Unit * fromUnit, Mantid::Kernel::Unit * toUnit); void convertUnitsViaTof(Mantid::Kernel::Unit * fromUnit, Mantid::Kernel::Unit * toUnit);
void convertUnitsQuickly(const double& factor, const double& power); void convertUnitsQuickly(const double& factor, const double& power);
template< class T >
void duplicateHelper(EventList* output, typename std::vector<T> & events) const;
void duplicate(EventList* output) const;
private: private:
///List of TofEvent (no weights). ///List of TofEvent (no weights).
......
...@@ -6,7 +6,6 @@ ...@@ -6,7 +6,6 @@
#include "MantidKernel/TimeSplitter.h" #include "MantidKernel/TimeSplitter.h"
#include "MantidDataObjects/TableWorkspace.h" #include "MantidDataObjects/TableWorkspace.h"
namespace Mantid namespace Mantid
{ {
namespace DataObjects namespace DataObjects
...@@ -57,6 +56,8 @@ namespace DataObjects ...@@ -57,6 +56,8 @@ namespace DataObjects
}; };
typedef boost::shared_ptr<SplittersWorkspace> SplittersWorkspace_sptr;
typedef boost::shared_ptr<const SplittersWorkspace> SplittersWorkspace_const_sptr;
} // namespace DataObjects } // namespace DataObjects
} // namespace Mantid } // namespace Mantid
......
...@@ -49,10 +49,12 @@ namespace DataObjects ...@@ -49,10 +49,12 @@ namespace DataObjects
} }
/** Compare two events' FRAME id, return true if e1 should be before e2. /** Compare two events' FRAME id, return true if e1 should be before e2.
* Assuming that if e1's pulse time is earlier than e2's, then e1 must be earlier regardless TOF value
* @param e1 :: first event * @param e1 :: first event
* @param e2 :: second event * @param e2 :: second event
* */ * */
bool compareEventPulseTimeTOF(const TofEvent& e1, const TofEvent& e2){ bool compareEventPulseTimeTOF(const TofEvent& e1, const TofEvent& e2)
{
if (e1.pulseTime() < e2.pulseTime()){ if (e1.pulseTime() < e2.pulseTime()){
return true; return true;
...@@ -65,8 +67,6 @@ namespace DataObjects ...@@ -65,8 +67,6 @@ namespace DataObjects
} }
//========================================================================== //==========================================================================
// ---------------------- EventList stuff ---------------------------------- // ---------------------- EventList stuff ----------------------------------
//========================================================================== //==========================================================================
...@@ -3262,11 +3262,10 @@ namespace DataObjects ...@@ -3262,11 +3262,10 @@ namespace DataObjects
} }
//------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------
/** Split the event list into n outputs, operating on a vector of either TofEvent's or WeightedEvent's /** Split the event list into n outputs, operating on a vector of either TofEvent's or WeightedEvent's
* Only event's pulse time is used to compare with splitters.
* It is a faster and simple version of splitByFullTimeHelper
* *
* @param splitter :: a TimeSplitterType giving where to split * @param splitter :: a TimeSplitterType giving where to split
* @param outputs :: a vector of where the split events will end up. The # of entries in there should * @param outputs :: a vector of where the split events will end up. The # of entries in there should
...@@ -3327,7 +3326,6 @@ namespace DataObjects ...@@ -3327,7 +3326,6 @@ namespace DataObjects
//Done! //Done!
} }
//------------------------------------------------------------------------------------------------ //------------------------------------------------------------------------------------------------
/** Split the event list into n outputs /** Split the event list into n outputs
* *
...@@ -3373,6 +3371,191 @@ namespace DataObjects ...@@ -3373,6 +3371,191 @@ namespace DataObjects
} }
//------------------------------------------------------------------------------------------------
/** Split the event list into n outputs, operating on a vector of either TofEvent's or WeightedEvent's
* The comparison between neutron event and splitter is based on neutron event's pulse time plus
*
* @param splitter :: a TimeSplitterType giving where to split
* @param outputs :: a vector of where the split events will end up. The # of entries in there should
* be big enough to accommodate the indices.
* @param events :: either this->events or this->weightedEvents.
*/
template< class T >
void EventList::splitByFullTimeHelper(Kernel::TimeSplitterType & splitter, std::map<int, EventList * > outputs,
typename std::vector<T> & events, double tofcorrection) const
{
// 1. Prepare to Iterate through the splitter at the same time
Kernel::TimeSplitterType::iterator itspl = splitter.begin();
Kernel::TimeSplitterType::iterator itspl_end = splitter.end();
int64_t start, stop;
int index;
// 2. Prepare to Iterate through all events (sorted by tof)
typename std::vector<T>::iterator itev = events.begin();
typename std::vector<T>::iterator itev_end = events.end();
// 3. This is the time of the first section. Anything before is thrown out.
while (itspl != itspl_end)
{
//Get the splitting interval times and destination
start = itspl->start().totalNanoseconds();
stop = itspl->stop().totalNanoseconds();
index = itspl->index();
// a) Skip the events before the start of the time
// TODO This step can be
EventList* myOutput = outputs[-1];
while (itev != itev_end)
{
int64_t fulltime = itev->m_pulsetime.totalNanoseconds() + static_cast<int64_t>(itev->m_tof*1000*tofcorrection);
if (fulltime < start)
{
// a1) Record to index = -1 space
const T eventCopy(*itev);
myOutput->addEventQuickly(eventCopy);
itev++;
}
else
{
break;
}
}
// b) Go through all the events that are in the interval (if any)
while (itev != itev_end)
{
int64_t fulltime = itev->m_pulsetime.totalNanoseconds() + static_cast<int64_t>(itev->m_tof*1000*tofcorrection);
if (fulltime < stop)
{
// b1) Copy the event into another
const T eventCopy(*itev);
EventList * myOutput = outputs[index];
//Add the copy to the output
myOutput->addEventQuickly(eventCopy);
++itev;
}
else
{
break;
}
}
//Go to the next interval
itspl++;
//But if we reached the end, then we are done.
if (itspl==itspl_end)
break;
//No need to keep looping through the filter if we are out of events
if (itev == itev_end)
break;
} // END-WHILE Splitter
return;
}
//------------------------------------------------------------------------------------------------
/** Split the event list into n outputs by event's full time (tof + pulse time)
*
* @param splitter :: a TimeSplitterType giving where to split
* @param outputs :: a map of where the split events will end up. The # of entries in there should
* be big enough to accommodate the indices.
* @param tofcorrection: a correction for each TOF to multiply with.
*/
void EventList::splitByFullTime(Kernel::TimeSplitterType & splitter, std::map<int, EventList * > outputs, double tofcorrection) const
{
if (eventType == WEIGHTED_NOTIME)
throw std::runtime_error("EventList::splitByTime() called on an EventList that no longer has time information.");
// 1. Start by sorting the event list by pulse time.
this->sortPulseTimeTOF();
// 2. Initialize all the outputs
std::map<int, EventList* >::iterator outiter;
for (outiter = outputs.begin(); outiter != outputs.end(); ++outiter)
{
EventList* opeventlist = outiter->second;
opeventlist->clear();
opeventlist->detectorIDs = this->detectorIDs;
opeventlist->refX = this->refX;
// Match the output event type.
opeventlist->switchTo(eventType);
}
//Do nothing if there are no entries
if (splitter.size() <= 0)
{
// 3A. Copy all events to group workspace = -1
this->duplicate(outputs[-1]);
}
else
{
// 3B. Split
switch (eventType)
{
case TOF:
splitByFullTimeHelper(splitter, outputs, this->events, tofcorrection);
break;
case WEIGHTED:
splitByFullTimeHelper(splitter, outputs, this->weightedEvents, tofcorrection);
break;
case WEIGHTED_NOTIME:
break;
}
}
return;
}
/*
* Duplicate helper function
*/
template< class T >
void EventList::duplicateHelper(EventList* output, typename std::vector<T> & events) const
{
for (size_t ie = 0; ie < events.size(); ++ ie)
{
const T eventCopy(events[ie]);
//Add the copy to the output
output->addEventQuickly(eventCopy);
}
return;
}
/*
* Duplicate events
*/
void EventList::duplicate(EventList* output) const
{
// 1. Initialize all the outputs
output->clear();
output->detectorIDs = this->detectorIDs;
output->refX = this->refX;
// Match the output event type.
output->switchTo(eventType);
// 2. Call Helper to do duplicate
switch (eventType)
{
case TOF:
duplicateHelper(output, this->events);
break;
case WEIGHTED:
duplicateHelper(output, this->weightedEvents);
break;
case WEIGHTED_NOTIME:
duplicateHelper(output, this->weightedEventsNoTime);
break;
}
return;
}
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** Get the vector of events contained in an EventList; /** Get the vector of events contained in an EventList;
* this is overloaded by event type. * this is overloaded by event type.
...@@ -3415,9 +3598,6 @@ namespace DataObjects ...@@ -3415,9 +3598,6 @@ namespace DataObjects
events = &el.getWeightedEventsNoTime(); events = &el.getWeightedEventsNoTime();
} }
//-------------------------------------------------------------------------- //--------------------------------------------------------------------------
/** Helper function for the conversion to TOF. This handles the different /** Helper function for the conversion to TOF. This handles the different
* event types. * event types.
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#include "MantidKernel/System.h" #include "MantidKernel/System.h"
#include "MantidAPI/Column.h" #include "MantidAPI/Column.h"
#include "MantidAPI/TableRow.h" #include "MantidAPI/TableRow.h"
#include "MantidKernel/IPropertyManager.h"
using namespace Mantid::Kernel; using namespace Mantid::Kernel;
using namespace Mantid::API; using namespace Mantid::API;
...@@ -85,3 +86,45 @@ namespace DataObjects ...@@ -85,3 +86,45 @@ namespace DataObjects
} // namespace Mantid } // namespace Mantid
} // namespace DataObjects } // namespace DataObjects
///\cond TEMPLATE
namespace Mantid
{
namespace Kernel
{
template<> DLLExport
Mantid::DataObjects::SplittersWorkspace_sptr IPropertyManager::getValue<Mantid::DataObjects::SplittersWorkspace_sptr>(const std::string &name) const
{
PropertyWithValue<Mantid::DataObjects::SplittersWorkspace_sptr>* prop =
dynamic_cast<PropertyWithValue<Mantid::DataObjects::SplittersWorkspace_sptr>*>(getPointerToProperty(name));
if (prop)
{
return *prop;
}
else
{
std::string message = "Attempt to assign property "+ name +" to incorrect type. Expected SplittersWorkspace.";
throw std::runtime_error(message);
}
}
template<> DLLExport
Mantid::DataObjects::SplittersWorkspace_const_sptr IPropertyManager::getValue<Mantid::DataObjects::SplittersWorkspace_const_sptr>(const std::string &name) const
{
PropertyWithValue<Mantid::DataObjects::SplittersWorkspace_sptr>* prop =
dynamic_cast<PropertyWithValue<Mantid::DataObjects::SplittersWorkspace_sptr>*>(getPointerToProperty(name));
if (prop)
{
return prop->operator()();
}
else
{
std::string message = "Attempt to assign property "+ name +" to incorrect type. Expected const SplittersWorkspace.";
throw std::runtime_error(message);
}
}
} // namespace Kernel
} // namespace Mantid
...@@ -10,6 +10,7 @@ ...@@ -10,6 +10,7 @@
#include "MantidAPI/SpectraDetectorMap.h" #include "MantidAPI/SpectraDetectorMap.h"
#include "MantidAPI/MatrixWorkspace.h" #include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/WorkspaceGroup.h" #include "MantidAPI/WorkspaceGroup.h"
#include "MantidAPI/Algorithm.h"
namespace WorkspaceCreationHelper namespace WorkspaceCreationHelper
{ {
...@@ -113,6 +114,8 @@ namespace WorkspaceCreationHelper ...@@ -113,6 +114,8 @@ namespace WorkspaceCreationHelper
Mantid::API::MatrixWorkspace_sptr createProcessedInelasticWS(const std::vector<double> &L2, const std::vector<double> &ploar, const std::vector<double> &azimutal, Mantid::API::MatrixWorkspace_sptr createProcessedInelasticWS(const std::vector<double> &L2, const std::vector<double> &ploar, const std::vector<double> &azimutal,
size_t numBins=4,double Emin=-10,double Emax=10,double Ei=11); size_t numBins=4,double Emin=-10,double Emax=10,double Ei=11);
Mantid::DataObjects::EventWorkspace_sptr createEventWorkspace3(Mantid::DataObjects::EventWorkspace_const_sptr sourceWS, std::string wsname, Mantid::API::Algorithm* alg);
}; };
......
...@@ -13,6 +13,10 @@ ...@@ -13,6 +13,10 @@
#include "MantidKernel/TimeSeriesProperty.h" #include "MantidKernel/TimeSeriesProperty.h"
#include "MantidAPI/Run.h" #include "MantidAPI/Run.h"
#include "MantidKernel/UnitFactory.h"
#include "MantidAPI/IAlgorithm.h"
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/SpectraAxis.h" #include "MantidAPI/SpectraAxis.h"
#include "MantidAPI/NumericAxis.h" #include "MantidAPI/NumericAxis.h"
...@@ -814,4 +818,75 @@ namespace WorkspaceCreationHelper ...@@ -814,4 +818,75 @@ namespace WorkspaceCreationHelper
} }
/*
* Create an EventWorkspace from a source EventWorkspace.
* The new workspace should be exactly the same as the source workspace but without any events
*/
Mantid::DataObjects::EventWorkspace_sptr createEventWorkspace3(Mantid::DataObjects::EventWorkspace_const_sptr sourceWS, std::string wsname, API::Algorithm* alg)
{
// 1. Initialize:use dummy numbers for arguments, for event workspace it doesn't matter
Mantid::DataObjects::EventWorkspace_sptr outputWS = Mantid::DataObjects::EventWorkspace_sptr(new DataObjects::EventWorkspace());
outputWS->setName(wsname);
outputWS->initialize(1,1,1);
// 2. Set the units
outputWS->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
outputWS->setYUnit("Counts");
outputWS->setTitle("Empty_Title");
// 3. Add the run_start property:
int runnumber = sourceWS->getRunNumber();
outputWS->mutableRun().addProperty("run_number", runnumber);
std::string runstartstr = sourceWS->run().getProperty("run_start")->value();
outputWS->mutableRun().addProperty("run_start", runstartstr);
// 4. Instrument
Mantid::API::Algorithm_sptr loadInst= alg->createSubAlgorithm("LoadInstrument");
// Now execute the sub-algorithm. Catch and log any error, but don't stop.
loadInst->setPropertyValue("InstrumentName", sourceWS->getInstrument()->getName());
loadInst->setProperty<MatrixWorkspace_sptr> ("Workspace", outputWS);
loadInst->setProperty("RewriteSpectraMap", true);
loadInst->executeAsSubAlg();
// Populate the instrument parameters in this workspace - this works around a bug
outputWS->populateInstrumentParameters();
// 6. Build spectrum and event list
// a) We want to pad out empty pixels.
detid2det_map detector_map;
outputWS->getInstrument()->getDetectors(detector_map);
// b) determine maximum pixel id
detid2det_map::iterator it;
detid_t detid_max = 0; // seems like a safe lower bound
for (it = detector_map.begin(); it != detector_map.end(); ++it)
if (it->first > detid_max)
detid_max = it->first;
// c) Pad all the pixels and Set to zero
std::vector<std::size_t> pixel_to_wkspindex;
pixel_to_wkspindex.reserve(detid_max+1); //starting at zero up to and including detid_max
pixel_to_wkspindex.assign(detid_max+1, 0);
size_t workspaceIndex = 0;
for (it = detector_map.begin(); it != detector_map.end(); ++it)
{
if (!it->second->isMonitor())
{
pixel_to_wkspindex[it->first] = workspaceIndex;
DataObjects::EventList & spec = outputWS->getOrAddEventList(workspaceIndex);
spec.addDetectorID(it->first);
// Start the spectrum number at 1
spec.setSpectrumNo(specid_t(workspaceIndex+1));
workspaceIndex += 1;
}
}
outputWS->doneAddingEventLists();
// Clear
pixel_to_wkspindex.clear();
return outputWS;
}
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment