Skip to content

Commit

Permalink
Merge pull request #14018 from mantidproject/13952_Fix_period_subtrac…
Browse files Browse the repository at this point in the history 
…tion_in_MuonCalculateAsymmetry

Fix period subtraction in MuonCalculateAsymmetry
  • Loading branch information
@KarlPalmen
KarlPalmen committed Oct 20, 2015
2 parents 570d45b + fc1c5e3 commit 08a3350
Show file tree
Hide file tree
Showing 5 changed files with 446 additions and 94 deletions.
22 changes: 16 additions & 6 deletions Framework/WorkflowAlgorithms/inc/MantidWorkflowAlgorithms/MuonCalculateAsymmetry.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -49,12 +49,22 @@ class DLLExport MuonCalculateAsymmetry : public API::Algorithm {
void init();
void exec();

/// Converts given workspace according to the OutputType
API::MatrixWorkspace_sptr convertWorkspace(API::MatrixWorkspace_sptr ws);

/// Merges two period workspaces according to PeriodOperation specified
API::MatrixWorkspace_sptr mergePeriods(API::MatrixWorkspace_sptr ws1,
API::MatrixWorkspace_sptr ws2);
// Calculates raw counts
API::MatrixWorkspace_sptr
calculateGroupCounts(const API::MatrixWorkspace_sptr &firstPeriodWS,
const API::MatrixWorkspace_sptr &secondPeriodWS,
int groupIndex, std::string op);
// Calculates asymmetry for specified spectrum
API::MatrixWorkspace_sptr
calculateGroupAsymmetry(const API::MatrixWorkspace_sptr &firstPeriodWS,
const API::MatrixWorkspace_sptr &secondPeriodWS,
int groupIndex, std::string op);
// Calculates asymmetry for a pair of spectra
API::MatrixWorkspace_sptr
calculatePairAsymmetry(const API::MatrixWorkspace_sptr &firstPeriodWS,
const API::MatrixWorkspace_sptr &secondPeriodWS,
int firstPairIndex, int secondPairIndex, double alpha,
std::string op);
};

} // namespace WorkflowAlgorithms
Expand Down
327 changes: 244 additions & 83 deletions Framework/WorkflowAlgorithms/src/MuonCalculateAsymmetry.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,125 +107,286 @@ void MuonCalculateAsymmetry::exec() {
MatrixWorkspace_sptr firstPeriodWS = getProperty("FirstPeriodWorkspace");
MatrixWorkspace_sptr secondPeriodWS = getProperty("SecondPeriodWorkspace");

MatrixWorkspace_sptr convertedWS;
// The type of calculation
const std::string type = getPropertyValue("OutputType");

if (secondPeriodWS) {
// Two periods
MatrixWorkspace_sptr mergedWS = mergePeriods(firstPeriodWS, secondPeriodWS);
convertedWS = convertWorkspace(mergedWS);
// The group index
int groupIndex = getProperty("GroupIndex");

// The type of period operation (+ or -)
std::string op = getProperty("PeriodOperation");

if (type == "GroupCounts") {

auto outWS =
calculateGroupCounts(firstPeriodWS, secondPeriodWS, groupIndex, op);

setProperty("OutputWorkspace", outWS);

} else if (type == "GroupAsymmetry") {

auto outWS =
calculateGroupAsymmetry(firstPeriodWS, secondPeriodWS, groupIndex, op);

setProperty("OutputWorkspace", outWS);

} else if (type == "PairAsymmetry") {

int pairFirstIndex = getProperty("PairFirstIndex");
int pairSecondIndex = getProperty("PairSecondIndex");
double alpha = getProperty("Alpha");

auto outWS =
calculatePairAsymmetry(firstPeriodWS, secondPeriodWS, pairFirstIndex,
pairSecondIndex, alpha, op);

setProperty("OutputWorkspace", outWS);

} else {
// Single period only
convertedWS = convertWorkspace(firstPeriodWS);
}

setProperty("OutputWorkspace", convertedWS);
throw std::invalid_argument("Specified OutputType is not supported");
}
}

/**
* Converts given workspace according to the OutputType.
* @param ws :: Workspace to convert
* @return Converted workspace
*/
MatrixWorkspace_sptr
MuonCalculateAsymmetry::convertWorkspace(MatrixWorkspace_sptr ws) {
const std::string type = getPropertyValue("OutputType");
* Calculates raw counts according to period operation
* @param firstPeriodWS :: [input] First period workspace
* @param secondPeriodWS :: [input] Second period workspace
* @param groupIndex :: [input] Index of the workspace to extract counts from
* @param op :: [input] Period operation (+ or -)
*/
MatrixWorkspace_sptr MuonCalculateAsymmetry::calculateGroupCounts(
const MatrixWorkspace_sptr &firstPeriodWS,
const MatrixWorkspace_sptr &secondPeriodWS, int groupIndex,
std::string op) {

if (type == "GroupCounts" || type == "GroupAsymmetry") {
int groupIndex = getProperty("GroupIndex");
if (secondPeriodWS) {
// Two periods supplied

if (groupIndex == EMPTY_INT())
throw std::runtime_error("GroupIndex is not specified");
MatrixWorkspace_sptr tempWS;

// Yank out the counts of requested group
IAlgorithm_sptr alg = createChildAlgorithm("ExtractSingleSpectrum");
alg->initialize();
alg->setProperty("InputWorkspace", ws);
alg->setProperty("WorkspaceIndex", groupIndex);
alg->execute();
if (op == "+") {

MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");
IAlgorithm_sptr alg = createChildAlgorithm("Plus");
alg->initialize();
alg->setProperty("LHSWorkspace", firstPeriodWS);
alg->setProperty("RHSWorkspace", secondPeriodWS);
alg->execute();
tempWS = alg->getProperty("OutputWorkspace");

if (type == "GroupAsymmetry") {
// GroupAsymmetry - counts with ExpDecay removed and normalized
} else if (op == "-") {

IAlgorithm_sptr alg = createChildAlgorithm("RemoveExpDecay");
IAlgorithm_sptr alg = createChildAlgorithm("Minus");
alg->initialize();
alg->setProperty("InputWorkspace", outWS);
alg->setProperty("LHSWorkspace", firstPeriodWS);
alg->setProperty("RHSWorkspace", secondPeriodWS);
alg->execute();

outWS = alg->getProperty("OutputWorkspace");
tempWS = alg->getProperty("OutputWorkspace");
}

return outWS;
} else if (type == "PairAsymmetry") {
// PairAsymmetry - result of AsymmetryCalc algorithm

int pairFirstIndex = getProperty("PairFirstIndex");
int pairSecondIndex = getProperty("PairSecondIndex");

if (pairFirstIndex == EMPTY_INT() || pairSecondIndex == EMPTY_INT())
throw std::invalid_argument("Both pair indices should be specified");

double alpha = getProperty("Alpha");
IAlgorithm_sptr alg = createChildAlgorithm("ExtractSingleSpectrum");
alg->initialize();
alg->setProperty("InputWorkspace", tempWS);
alg->setProperty("WorkspaceIndex", groupIndex);
alg->execute();
MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");

// We get pair groups as their workspace indices, but AsymmetryCalc wants
// spectra numbers,
// so need to convert
specid_t spectraNo1 = ws->getSpectrum(pairFirstIndex)->getSpectrumNo();
specid_t spectraNo2 = ws->getSpectrum(pairSecondIndex)->getSpectrumNo();
return outWS;

if (spectraNo1 == -1 || spectraNo2 == -1 || spectraNo1 == spectraNo2)
throw std::invalid_argument(
"Spectra numbers of the input workspace are not set properly");
} else {
// Only one period supplied

IAlgorithm_sptr alg = createChildAlgorithm("AsymmetryCalc");
alg->setProperty("InputWorkspace", ws);
// As strings, cause otherwise would need to create arrays with single
// elements
alg->setPropertyValue("ForwardSpectra",
boost::lexical_cast<std::string>(spectraNo1));
alg->setPropertyValue("BackwardSpectra",
boost::lexical_cast<std::string>(spectraNo2));
alg->setProperty("Alpha", alpha);
IAlgorithm_sptr alg = createChildAlgorithm("ExtractSingleSpectrum");
alg->initialize();
alg->setProperty("InputWorkspace", firstPeriodWS);
alg->setProperty("WorkspaceIndex", groupIndex);
alg->execute();

MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");

return outWS;
}

throw std::invalid_argument("Specified OutputType is not supported");
}

/**
* Merges two period workspaces according to PeriodOperation specified.
* @param ws1 :: First period workspace
* @param ws2 :: Second period workspace
* @return Merged workspace
*/
MatrixWorkspace_sptr
MuonCalculateAsymmetry::mergePeriods(MatrixWorkspace_sptr ws1,
MatrixWorkspace_sptr ws2) {
std::string op = getProperty("PeriodOperation");
* Calculates single-spectrum asymmetry according to period operation
* @param firstPeriodWS :: [input] First period workspace
* @param secondPeriodWS :: [input] Second period workspace
* @param groupIndex :: [input] Workspace index for which to calculate asymmetry
* @param op :: [input] Period operation (+ or -)
*/
MatrixWorkspace_sptr MuonCalculateAsymmetry::calculateGroupAsymmetry(
const MatrixWorkspace_sptr &firstPeriodWS,
const MatrixWorkspace_sptr &secondPeriodWS, int groupIndex,
std::string op) {

// The output workspace
MatrixWorkspace_sptr tempWS;

if (secondPeriodWS) {
// Two period workspaces where supplied

if (op == "+") {

std::string algorithmName;
// Sum
IAlgorithm_sptr alg = createChildAlgorithm("Plus");
alg->initialize();
alg->setProperty("LHSWorkspace", firstPeriodWS);
alg->setProperty("RHSWorkspace", secondPeriodWS);
alg->execute();
MatrixWorkspace_sptr sumWS = alg->getProperty("OutputWorkspace");

// GroupIndex as vector
// Necessary if we want RemoveExpDecay to fit only the requested
// spectrum
std::vector<int> spec(1, groupIndex);

// Remove decay
IAlgorithm_sptr asym = createChildAlgorithm("RemoveExpDecay");
asym->setProperty("InputWorkspace", sumWS);
asym->setProperty("Spectra", spec);
asym->execute();
tempWS = asym->getProperty("OutputWorkspace");

} else if (op == "-") {

// GroupIndex as vector
std::vector<int> spec(1, groupIndex);

// Remove decay (first period ws)
IAlgorithm_sptr asym = createChildAlgorithm("RemoveExpDecay");
asym->initialize();
asym->setProperty("InputWorkspace", firstPeriodWS);
asym->setProperty("Spectra", spec);
asym->execute();
MatrixWorkspace_sptr asymFirstPeriod =
asym->getProperty("OutputWorkspace");

// Remove decay (second period ws)
asym = createChildAlgorithm("RemoveExpDecay");
asym->initialize();
asym->setProperty("InputWorkspace", secondPeriodWS);
asym->setProperty("Spectra", spec);
asym->execute();
MatrixWorkspace_sptr asymSecondPeriod =
asym->getProperty("OutputWorkspace");

// Now subtract
IAlgorithm_sptr alg = createChildAlgorithm("Minus");
alg->initialize();
alg->setProperty("LHSWorkspace", asymFirstPeriod);
alg->setProperty("RHSWorkspace", asymSecondPeriod);
alg->execute();
tempWS = alg->getProperty("OutputWorkspace");
}
} else {
// Only one period was supplied

if (op == "+") {
algorithmName = "Plus";
} else if (op == "-") {
algorithmName = "Minus";
IAlgorithm_sptr alg = createChildAlgorithm("RemoveExpDecay");
alg->initialize();
alg->setProperty("InputWorkspace", firstPeriodWS);
alg->execute();
tempWS = alg->getProperty("OutputWorkspace");
}

IAlgorithm_sptr alg = createChildAlgorithm(algorithmName);
// Extract the requested spectrum
IAlgorithm_sptr alg = createChildAlgorithm("ExtractSingleSpectrum");
alg->initialize();
alg->setProperty("LHSWorkspace", ws1);
alg->setProperty("RHSWorkspace", ws2);
alg->setProperty("InputWorkspace", tempWS);
alg->setProperty("WorkspaceIndex", groupIndex);
alg->execute();

MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");

return outWS;
}

/**
* Calculates pair asymmetry according to period operation
* @param firstPeriodWS :: [input] First period workspace
* @param secondPeriodWS :: [input] Second period workspace
* @param firstPairIndex :: [input] Workspace index for the forward group
* @param secondPairIndex :: [input] Workspace index for the backward group
* @param alpha :: [input] The balance parameter
* @param op :: [input] Period operation (+ or -)
*/
MatrixWorkspace_sptr MuonCalculateAsymmetry::calculatePairAsymmetry(
const MatrixWorkspace_sptr &firstPeriodWS,
const MatrixWorkspace_sptr &secondPeriodWS, int firstPairIndex,
int secondPairIndex, double alpha, std::string op) {

// Pair indices as vectors
std::vector<int> fwd(1, firstPairIndex + 1);
std::vector<int> bwd(1, secondPairIndex + 1);

if (secondPeriodWS) {

MatrixWorkspace_sptr outWS;

if (op == "+") {

// Sum
IAlgorithm_sptr alg = createChildAlgorithm("Plus");
alg->initialize();
alg->setProperty("LHSWorkspace", firstPeriodWS);
alg->setProperty("RHSWorkspace", secondPeriodWS);
alg->execute();
MatrixWorkspace_sptr sumWS = alg->getProperty("OutputWorkspace");

// Asymmetry calculation
alg = createChildAlgorithm("AsymmetryCalc");
alg->setProperty("InputWorkspace", sumWS);
alg->setProperty("ForwardSpectra", fwd);
alg->setProperty("BackwardSpectra", bwd);
alg->setProperty("Alpha", alpha);
alg->execute();
outWS = alg->getProperty("OutputWorkspace");

} else if (op == "-") {

std::vector<int> fwd(1, firstPairIndex + 1);
std::vector<int> bwd(1, secondPairIndex + 1);

// First period asymmetry
IAlgorithm_sptr alg = createChildAlgorithm("AsymmetryCalc");
alg->setProperty("InputWorkspace", firstPeriodWS);
alg->setProperty("ForwardSpectra", fwd);
alg->setProperty("BackwardSpectra", bwd);
alg->setProperty("Alpha", alpha);
alg->execute();
MatrixWorkspace_sptr asymFirstPeriod =
alg->getProperty("OutputWorkspace");

// Second period asymmetry
alg = createChildAlgorithm("AsymmetryCalc");
alg->setProperty("InputWorkspace", secondPeriodWS);
alg->setProperty("ForwardSpectra", fwd);
alg->setProperty("BackwardSpectra", bwd);
alg->setProperty("Alpha", alpha);
alg->execute();
MatrixWorkspace_sptr asymSecondPeriod =
alg->getProperty("OutputWorkspace");

// Now subtract
alg = createChildAlgorithm("Minus");
alg->initialize();
alg->setProperty("LHSWorkspace", asymFirstPeriod);
alg->setProperty("RHSWorkspace", asymSecondPeriod);
alg->execute();
outWS = alg->getProperty("OutputWorkspace");
}

return outWS;

} else {

IAlgorithm_sptr alg = createChildAlgorithm("AsymmetryCalc");
alg->setProperty("InputWorkspace", firstPeriodWS);
alg->setProperty("ForwardSpectra", fwd);
alg->setProperty("BackwardSpectra", bwd);
alg->setProperty("Alpha", alpha);
alg->execute();
MatrixWorkspace_sptr outWS = alg->getProperty("OutputWorkspace");

return outWS;
}
}
} // namespace WorkflowAlgorithms
} // namespace Mantid
Loading

0 comments on commit 08a3350

Please sign in to comment.