Unverified Commit 746b0260 authored by Zhou, Wenduo's avatar Zhou, Wenduo Committed by GitHub
Browse files

Merge pull request #31541 from mantidproject/SPEC94_detailed_balance_31459_clean

PR is merged as https://github.com/mantidproject/mantid/pull/31536 has been merged.
parents b14cdd12 57b7ee90
......@@ -2,6 +2,7 @@
set(SRC_FILES
src/AccumulateMD.cpp
src/AndMD.cpp
src/ApplyDetailedBalanceMD.cpp
src/BaseConvertToDiffractionMDWorkspace.cpp
src/BinMD.cpp
src/BinaryOperationMD.cpp
......@@ -117,6 +118,7 @@ set(
INC_FILES
inc/MantidMDAlgorithms/AccumulateMD.h
inc/MantidMDAlgorithms/AndMD.h
inc/MantidMDAlgorithms/ApplyDetailedBalanceMD.h
inc/MantidMDAlgorithms/BaseConvertToDiffractionMDWorkspace.h
inc/MantidMDAlgorithms/BinMD.h
inc/MantidMDAlgorithms/BinaryOperationMD.h
......@@ -234,6 +236,7 @@ set(TEST_FILES
# these tests are as they test verify different parts of the CPR algorithms
AccumulateMDTest.h
AndMDTest.h
ApplyDetailedBalanceMDTest.h
BooleanBinaryOperationMDTest.h
CalculateCoverageDGSTest.h
CentroidPeaksMD2Test.h
......
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2011 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once
#include "MantidAPI/Algorithm.h"
#include "MantidAPI/IMDEventWorkspace_fwd.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include "MantidMDAlgorithms/DllConfig.h"
namespace Mantid {
namespace MDAlgorithms {
/** ApplyDetailedBalance : Perform the And boolean operation on two MDHistoWorkspaces
@date 2011-11-08
*/
class MANTID_MDALGORITHMS_DLL ApplyDetailedBalanceMD : public API::Algorithm {
public:
ApplyDetailedBalanceMD() : mDeltaEIndex(999) {}
const std::string name() const override;
int version() const override;
const std::string category() const override;
const std::string summary() const override;
private:
/// Initialize the proeprties
void init() override;
/// Run the algorithm
void exec() override;
/// Validate inputs
std::map<std::string, std::string> validateInputs() override;
/// Apply detailed balance to each MDEvent
template <typename MDE, size_t nd>
void applyDetailedBalance(typename Mantid::DataObjects::MDEventWorkspace<MDE, nd>::sptr ws);
/// Get temperature
std::string getTemperature(API::IMDEventWorkspace_sptr mdws);
/// Check input workspace dimension
std::string checkInputMDDimension();
/// index of the MD dimension index for DeltaE
size_t mDeltaEIndex;
/// map of temperature retrieved from sample logs
std::map<uint16_t, double> mExpinfoTemperatureMean;
};
} // namespace MDAlgorithms
} // namespace Mantid
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidMDAlgorithms/ApplyDetailedBalanceMD.h"
#include "MantidDataObjects/MDBoxIterator.h"
#include "MantidDataObjects/MDEventFactory.h"
// #include "MantidAPI/IMDEventWorkspace.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/MDGeometry.h"
#include "MantidAPI/MultipleExperimentInfos.h"
#include "MantidAPI/Run.h"
#include "MantidGeometry/MDGeometry/IMDDimension.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/Property.h"
#include "MantidKernel/PropertyWithValue.h"
#include "MantidKernel/SpecialCoordinateSystem.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidMDAlgorithms/MDWSDescription.h"
#include <math.h>
using namespace Mantid::Kernel;
using namespace Mantid::API;
using namespace Mantid::Geometry;
using namespace Mantid::DataObjects;
namespace Mantid {
namespace MDAlgorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(ApplyDetailedBalanceMD)
/// Algorithm's name for identification. @see Algorithm::name
const std::string ApplyDetailedBalanceMD::name() const { return "ApplyDetailedBalanceMD"; }
/// Algorithm's version for identification. @see Algorithm::version
int ApplyDetailedBalanceMD::version() const { return 1; }
/// Summary
const std::string ApplyDetailedBalanceMD::summary() const { return "Apply detailed balance to MDEventWorkspace"; }
/// category
const std::string ApplyDetailedBalanceMD::category() const { return "MDAlgorithms"; }
//---------------------------------------------------------------------------------------------------------
/**
* @brief Define input and output properties
*/
void ApplyDetailedBalanceMD::init() {
declareProperty(
std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("InputWorkspace", "", Kernel::Direction::Input),
"An input MDEventWorkspace. Must be in Q_sample/Q_lab frame. Must have an axis as DeltaE");
declareProperty(std::make_unique<PropertyWithValue<std::string>>("Temperature", "", Direction::Input),
"SampleLog variable name that contains the temperature or a number");
declareProperty(
std::make_unique<WorkspaceProperty<API::IMDEventWorkspace>>("OutputWorkspace", "", Kernel::Direction::Output),
"The output MDEventWorkspace with detailed balance applied");
}
//---------------------------------------------------------------------------------------------------------
/**
* @brief Main execution body
*/
void ApplyDetailedBalanceMD::exec() {
// Get input workspace
API::IMDEventWorkspace_sptr input_ws = getProperty("InputWorkspace");
// Process input workspace and create output workspace
std::string output_ws_name = getPropertyValue("OutputWorkspace");
API::IMDEventWorkspace_sptr output_ws(0);
if (input_ws->getName() == output_ws_name) {
// Calcualte in-place
output_ws = input_ws;
} else {
// Clone input workace to output workspace
output_ws = input_ws->clone();
}
// Apply detailed balance to MDEvents
CALL_MDEVENT_FUNCTION(applyDetailedBalance, output_ws);
// refresh cache for MDBoxes: set correct Box signal
output_ws->refreshCache();
// Clear masking (box flags) from the output workspace
output_ws->clearMDMasking();
// Set output
setProperty("OutputWorkspace", output_ws);
}
//---------------------------------------------------------------------------------------------------------
/**
* @brief Validate inputs
* Input MDEventWorkspace dimensions:
* - in Q_sample or Q_lab frame, the 4th dimension is DeltaE, or
* - the first dimension is |Q| and second is DeltaE
* validate that run objects have Ei defined (number greater than 0)
* first input is an MD event workspace as validated above
second input is Temperature - must be either a float>0, or a string pointing to a log name.
if temperature points to a log name, it must be present in each experiment info, and the average must be greater than 0
*/
std::map<std::string, std::string> ApplyDetailedBalanceMD::validateInputs() {
std::map<std::string, std::string> output;
// Get input workspace
API::IMDEventWorkspace_sptr input_ws = getProperty("InputWorkspace");
// check input dimension
std::string dim_error = checkInputMDDimension();
if (dim_error != "") {
output["InputWorkspace"] = dim_error;
}
std::string kerror = getTemperature(input_ws);
if (kerror != "")
output["Temperature"] = kerror;
return output;
}
//---------------------------------------------------------------------------------------------------------
/**
* @brief Check input MDEventWorkspace dimension
* validate dimensions: input workspace is in Q_sample or Q_lab frame, and the 4th dimension is DeltaE, or the first
* dimension is |Q| and second is DeltaE
* @return
*/
std::string ApplyDetailedBalanceMD::checkInputMDDimension() {
std::string errormsg("");
API::IMDEventWorkspace_sptr inputws = getProperty("InputWorkspace");
size_t numdims = inputws->getNumDims();
// Get and check the dimensions: Q3D or Q1D
const Mantid::Kernel::SpecialCoordinateSystem coordsys = inputws->getSpecialCoordinateSystem();
size_t qdim(0);
std::string qdimstr("Not Q3D or |Q|");
if (coordsys == Mantid::Kernel::SpecialCoordinateSystem::QLab ||
coordsys == Mantid::Kernel::SpecialCoordinateSystem::QSample) {
// q3d
qdim = 3;
qdimstr = "Q3D";
} else {
// search Q1D: at any place
for (size_t i = 0; i < numdims; ++i) {
if (inputws->getDimension(i)->getName() == "|Q|") {
qdim = 1;
qdimstr = "|Q|";
break;
}
}
}
// Check DeltaE
if (qdim == 1 && inputws->getDimension(1)->getName() == "DeltaE") {
// 2nd dimension is DeltaE
mDeltaEIndex = 1;
} else if (qdim == 3 && inputws->getDimension(3)->getName() == "DeltaE") {
// 4th dimension is DeltaE
mDeltaEIndex = 3;
} else {
// Error
g_log.error() << "Coordiate system = " << coordsys << " does not meet requirement: \n";
for (size_t i = 0; i < numdims; ++i) {
g_log.error() << i << "-th dim: " << inputws->getDimension(i)->getName() << "\n";
}
errormsg += "Q Dimension (" + qdimstr +
") is neither Q3D nor |Q|. Or DeltaE is found in proper place (2nd or 4th dimension).";
}
return errormsg;
}
//---------------------------------------------------------------------------------------------
/**
* @brief Apply detailed balance to each MDEvent in MDEventWorkspace
*/
template <typename MDE, size_t nd>
void ApplyDetailedBalanceMD::applyDetailedBalance(typename Mantid::DataObjects::MDEventWorkspace<MDE, nd>::sptr ws) {
// Get Box from MDEventWorkspace
MDBoxBase<MDE, nd> *box1 = ws->getBox();
std::vector<API::IMDNode *> boxes;
box1->getBoxes(boxes, 1000, true);
auto numBoxes = int(boxes.size());
// Add the boxes in parallel. They should be spread out enough on each
// core to avoid stepping on each other.
// cppcheck-suppress syntaxError
PRAGMA_OMP( parallel for if (!ws->isFileBacked()))
for (int i = 0; i < numBoxes; ++i) {
PARALLEL_START_INTERUPT_REGION
auto *box = dynamic_cast<MDBox<MDE, nd> *>(boxes[i]);
if (box && !box->getIsMasked()) {
// get the MEEvents from box
std::vector<MDE> &events = box->getEvents();
// Add events, with bounds checking
for (auto it = events.begin(); it != events.end(); ++it) {
// Create the event
// do calculattion
float temperatue(static_cast<float>(mExpinfoTemperatureMean[it->getExpInfoIndex()]));
// delta_e = it->getCenter(mDeltaEIndex);
// factor = pi * (1 - exp(-deltaE/(kb*T)))
float factor = static_cast<float>(M_PI) *
(static_cast<float>(1.) - exp(-it->getCenter(mDeltaEIndex) *
static_cast<float>(PhysicalConstants::meVtoKelvin / temperatue)));
// calcalate and set intesity
auto intensity = it->getSignal() * factor;
it->setSignal(intensity);
// calculate and set error
auto error2 = it->getErrorSquared() * factor * factor;
// error2 *= factor * factor;
it->setErrorSquared(error2);
}
}
box->releaseEvents();
PARALLEL_END_INTERUPT_REGION
}
PARALLEL_CHECK_INTERUPT_REGION
return;
}
//---------------------------------------------------------------------------------------------------------
/**
* @brief Retrieve sample temperature
* Temperature value can be specified by either property Temperature, or
* it can be calcualted from sample temperture log in the MDWorkspace
*/
std::string ApplyDetailedBalanceMD::getTemperature(API::IMDEventWorkspace_sptr mdws) {
// Get temperture sample log name
std::string Tstring = getProperty("Temperature");
std::string temperature_error("");
// Try to convert Tstring to a float
float temperature;
try {
temperature = boost::lexical_cast<float>(Tstring);
} catch (...) {
// set to a unphysical value
temperature = -10;
}
// the input property could be a valid float; if not must search the experiment info
mExpinfoTemperatureMean.clear();
uint16_t numexpinfo = mdws->getNumExperimentInfo();
for (uint16_t i = 0; i < numexpinfo; ++i) {
if (temperature < 0) {
// if user specified is not a valid float
ExperimentInfo_const_sptr expinfo = mdws->getExperimentInfo(i);
if (expinfo->run().hasProperty(Tstring)) {
auto log = dynamic_cast<Kernel::TimeSeriesProperty<double> *>(expinfo->run().getProperty(Tstring));
if (!log) {
// wrong type of sample log: must be TimeSeriesProperty<double>
std::stringstream errss;
errss << "ExperimentInfo" << i << " has " << Tstring << ", which is not a valid double-valuesd log";
temperature_error += errss.str() + "\n";
}
mExpinfoTemperatureMean[i] = log->getStatistics().mean;
} else {
// specified sample log does not exist
std::stringstream errss;
errss << "ExperimentInfo " << i << " does not have tempertaure log " << Tstring;
temperature_error += errss.str() + "\n";
}
} else {
// set user specified temperature to map
mExpinfoTemperatureMean[i] = temperature;
}
}
return temperature_error;
}
} // namespace MDAlgorithms
} // namespace Mantid
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright &copy; 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#pragma once
#include "MantidKernel/System.h"
#include "MantidKernel/Timer.h"
#include <cxxtest/TestSuite.h>
#include "MantidAPI/AlgorithmManager.h"
#include "MantidAPI/AnalysisDataService.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidDataHandling/MoveInstrumentComponent.h"
#include "MantidDataObjects/MDEventWorkspace.h"
#include "MantidDataObjects/MDHistoWorkspace.h"
#include "MantidKernel/PhysicalConstants.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidMDAlgorithms/ApplyDetailedBalanceMD.h"
#include "MantidMDAlgorithms/CompareMDWorkspaces.h"
#include "MantidMDAlgorithms/ConvertToMD.h"
#include "MantidMDAlgorithms/MergeMD.h"
using namespace Mantid;
using namespace Mantid::MDAlgorithms;
using namespace Mantid::API;
using Mantid::DataObjects::MDHistoWorkspace_sptr;
class ApplyDetailedBalanceMDTest : public CxxTest::TestSuite {
public:
void test_Init() {
ApplyDetailedBalanceMD alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize())
TS_ASSERT(alg.isInitialized())
}
//---------------------------------------------------------------------------------------------
/**
* @brief Test apply detailed balance to 1 run
*/
void test_1run() {
// Check whether the MD to test does exist
auto singlemd = Mantid::API::AnalysisDataService::Instance().retrieve(mMDWorkspace1Name);
TS_ASSERT(singlemd);
// specify the output
std::string outputname("DetailedBalanceSingleQ3Test");
// Calculate detailed balance for the single MDEventWorkspace
ApplyDetailedBalanceMD alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize());
TS_ASSERT(alg.isInitialized());
alg.setPropertyValue("InputWorkspace", mMDWorkspace1Name);
alg.setProperty("Temperature", "SampleTemp");
alg.setProperty("OutputWorkspace", outputname);
alg.execute();
TS_ASSERT(alg.isExecuted());
// Verify
TS_ASSERT(AnalysisDataService::Instance().doesExist(outputname));
bool equals = compareMDEvents(outputname, GoldDetailedBalanceSingleWSName);
TS_ASSERT(equals);
// Clean up
cleanWorkspace(outputname, false);
}
/**
* @brief Test applying detailed balance to 2 merged runs
*/
void test_merged_runs() {
auto mergedmd = Mantid::API::AnalysisDataService::Instance().retrieve(mMergedWorkspaceName);
TS_ASSERT(mergedmd);
// specify the output
std::string outputname("DetailedBalanceMergedQ3Test");
ApplyDetailedBalanceMD alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize());
TS_ASSERT(alg.isInitialized());
alg.setPropertyValue("InputWorkspace", mMergedWorkspaceName);
alg.setProperty("Temperature", "SampleTemp");
alg.setProperty("OutputWorkspace", outputname);
alg.execute();
TS_ASSERT(alg.isExecuted());
// Verify
TS_ASSERT(AnalysisDataService::Instance().doesExist(outputname));
bool equals = compareMDEvents(outputname, gold_detail_balanced_merged_name);
TS_ASSERT(equals);
// clean up
cleanWorkspace(outputname, false);
}
/**
* @brief Test input workspace with |Q| and without temperature
*/
void test_q1d_run() {
auto q1dmd = Mantid::API::AnalysisDataService::Instance().retrieve(mMDWorkspaceQ1Dname);
TS_ASSERT(q1dmd);
ApplyDetailedBalanceMD alg;
TS_ASSERT_THROWS_NOTHING(alg.initialize());
TS_ASSERT(alg.isInitialized());
alg.setPropertyValue("InputWorkspace", mMDWorkspaceQ1Dname);
alg.setProperty("Temperature", "SampleTemp");
alg.setProperty("OutputWorkspace", "OutputDetaiedBalanceQ1D");
// Expect to fail due to missing temperature sample log
TS_ASSERT_THROWS_ANYTHING(alg.execute());
// Set temperature explicitly
alg.setProperty("Temperature", "1.2345");
TS_ASSERT_THROWS_NOTHING(alg.execute());
// Check existence and clean
cleanWorkspace("OutputDetaiedBalanceQ1D", true);
}
/**
* @brief Clean workspace from ADS if it does exist
*/
void cleanWorkspace(const std::string &wsname, bool assert_existence) {
bool ws_exist = Mantid::API::AnalysisDataService::Instance().doesExist(wsname);
// assert existence
if (assert_existence)
TS_ASSERT(ws_exist);
// clean from ADS
if (ws_exist)
Mantid::API::AnalysisDataService::Instance().remove(wsname);
}
void setUp() override {
// Define workspace names
mEventWSName = "DetailedBalanceRawEvent";
mMDWorkspace1Name = "DetailedBalanceInputSingleMDEvent";
mMergedWorkspaceName = "DetailedBalanceInputMergedMDEvent";
mMDWorkspaceQ1Dname = "DetaledBalanceInputQ1DMDEvent";
// Gold workspace names
gold_detail_balanced_merged_name = "DetaildBalancedMergedGoldMD";
GoldDetailedBalanceSingleWSName = "DetailedBalanceSingleGoldMD";
// Prepare (first) sample workspace
createSampleWorkspace(mEventWSName);
// add sample log Ei
addSampleLog(mEventWSName, "Ei", "20.", "Number");
// move bank 1
moveBank(mEventWSName, "bank1", 3, 3);
// move bank 2
moveBank(mEventWSName, "bank2", -3, -3);
// ad sample log for temperature
addSampleLog(mEventWSName, "SampleTemp", "25.0", "Number Series");
// set geoniometer
setGoniometer(mEventWSName, "Axis0", "0,0,1,0,1");
// convert to MD
convertToMD(mEventWSName, mMDWorkspace1Name, "Q3D");
// Prepare the 2nd MDEventWorkspace
const std::string event_ws_name2("DetailedBalance2WS");
const std::string md_ws_name2("DetailedBalance2MD");
createSampleWorkspace(event_ws_name2);
// add sample log Ei
addSampleLog(event_ws_name2, "Ei", "20.", "Number");
// move bank 1
moveBank(event_ws_name2, "bank1", 3, 3);
// move bank 2
moveBank(event_ws_name2, "bank2", -3, -3);
// ad sample log for temperature
addSampleLog(event_ws_name2, "SampleTemp", "250.0", "Number Series");
// set geoniometer
setGoniometer(event_ws_name2, "Axis0", "30,0,1,0,1");
// convert to MD
convertToMD(event_ws_name2, md_ws_name2, "Q3D");
// Prepare the 3rd MDEventWorkspace for |Q| and without sample temperature
const std::string event_ws_name3("DetailedBalance3WS");
createSampleWorkspace(event_ws_name3);
// add sample log Ei
addSampleLog(event_ws_name3, "Ei", "20.", "Number");
// move bank 1
moveBank(event_ws_name3, "bank1", 3, 3);
// move bank 2
moveBank(event_ws_name3, "bank2", -3, -3);
// set geoniometer
setGoniometer(event_ws_name3, "Axis0", "30,0,1,0,1");
// convert to MD
convertToMD(event_ws_name3, mMDWorkspaceQ1Dname, "|Q|");
// Merge 2 workspace
std::string workspaces(mMDWorkspace1Name + ", " + md_ws_name2);
MergeMD merge_alg;
merge_alg.initialize();
merge_alg.setProperty("InputWorkspaces", workspaces);
merge_alg.setProperty("OutputWorkspace", mMergedWorkspaceName);
merge_alg.execute();
// Calculate the expected result from existing algorithms
calculate_detailed_balance(mEventWSName, event_ws_name2, GoldDetailedBalanceSingleWSName,
gold_detail_balanced_merged_name);
// clean the temporary workspaces
Mantid::API::AnalysisDataService::Instance().remove(event_ws_name2);
Mantid::API::AnalysisDataService::Instance().remove(md_ws_name2);
// FIXME: clean DetailedBalance3WS
// FIXME: clean DetailedBalance2WS
}
void tearDown() override {
// clean up
bool exist = Mantid::API::AnalysisDataService::Instance().doesExist(mEventWSName);
TS_ASSERT(exist);
Mantid::API::AnalysisDataService::Instance().remove(mEventWSName);
// single MD workspace
bool single_exist = Mantid::API::AnalysisDataService::Instance().doesExist(mMDWorkspace1Name);
TS_ASSERT(single_exist);
Mantid::API::AnalysisDataService::Instance().remove(mMDWorkspace1Name);
// merged MD workspace
bool merge_exist = Mantid::API::AnalysisDataService::Instance().doesExist(mMergedWorkspaceName);
TS_ASSERT(merge_exist);
Mantid::API::AnalysisDataService::Instance().remove(mMergedWorkspaceName);
}
private:
std::string mEventWSName;
std::string mMDWorkspace1Name;