Skip to content
Snippets Groups Projects
Commit 54de80c6 authored by Federico Montesino Pouzols's avatar Federico Montesino Pouzols
Browse files

Merge pull request #13686 from mantidproject/11588_PhaseQuad_improvements 

PhaseQuadMuon improvements
parents 1f14ffd9 48a575b2
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PhaseQuadMuon.h
+ 9
− 66
View file @ 54de80c6
......@@ -37,18 +37,14 @@ Code Documentation is available at: <http://doxygen.mantidproject.org>
class DLLExport PhaseQuadMuon : public API::Algorithm {
public:
/// Default constructor
PhaseQuadMuon()
: API::Algorithm(), m_muLife(2.19703), m_bigNumber(1e10), m_tPulseOver(0),
m_pulseTail(182), m_poissonLim(30), m_pulseTwo(0.325), m_nHist(0),
m_nData(0), m_res(0.0), m_meanLag(0.0), m_tValid(0), m_isDouble(false),
m_tMin(0.0){};
PhaseQuadMuon() {};
/// Destructor
virtual ~PhaseQuadMuon(){};
/// Algorithm's name for identification overriding a virtual method
virtual const std::string name() const { return "PhaseQuad"; }
/// Summary of algorithm's purpose
virtual const std::string summary() const {
return "Calculate Muon squashograms from InputWorkspace and PhaseTable/PhaseList.";
return "Generates a quadrature phase signal.";
}
/// Algorithm's version for identification overriding a virtual method
......@@ -57,71 +53,18 @@ public:
virtual const std::string category() const { return "Muon"; }
private:
class HistData {
public:
bool detOK; // Detector is OK
double n0; // Detector n0
double alpha; // Detector efficiency
double phi; // Detector phase
};
/// Initialise the properties
void init();
/// Run the algorithm
void exec();
/// Convert X units from micro-secs to nano-secs and shift to start at t=0
void convertToNanoSecs(API::MatrixWorkspace_sptr inputWs);
/// Convert X units from nano-secs to micro-secs and shift back
void convertToMicroSecs(API::MatrixWorkspace_sptr inputWs);
/// Load the Phase Table
void loadPhaseTable(API::ITableWorkspace_sptr phaseTable,
API::ITableWorkspace_sptr deadTimeTable);
/// Load the Phase List
void loadPhaseList(const std::string &filename,
API::ITableWorkspace_sptr deadTimeTable);
/// Apply dead time correction
void deadTimeCorrection(API::MatrixWorkspace_sptr inputWs,
API::ITableWorkspace_sptr deadTimeTable,
API::MatrixWorkspace_sptr &tempWs);
/// Rescale detector efficiency to maximum value
void normaliseAlphas(std::vector<HistData> &m_histData);
/// Remove exponential decay from input histograms
void loseExponentialDecay(API::MatrixWorkspace_sptr tempWs);
/// Validate inputs
std::map<std::string, std::string> validateInputs();
/// Calculate the normalization constants
std::vector<double> getExponentialDecay(const API::MatrixWorkspace_sptr &ws);
/// Create squashograms
void squash(const API::MatrixWorkspace_sptr tempWs,
API::MatrixWorkspace_sptr outputWs);
/// Put back in exponential decay
void regainExponential(API::MatrixWorkspace_sptr outputWs);
/// Muon lifetime
double m_muLife;
/// Maximum counts expected
double m_bigNumber;
/// Pulse definitely finished by here (bin no)
int m_tPulseOver;
/// Number of bins to exclude after lag-time (ie pulse arrival time)
double m_pulseTail;
/// Poisson limit
double m_poissonLim;
/// Time (microsec) by which a well-def'd point in the first proton/pion pulse
/// leads its counterpart in the second
double m_pulseTwo;
/// Number of input histograms
int m_nHist;
/// Number of datapoints per histogram
int m_nData;
/// Time resolution
double m_res;
/// Mean of time-shifts
double m_meanLag;
/// Good muons from here on (bin no). Unused now but can be needed in the
/// future
int m_tValid;
/// Double-pulse flag
bool m_isDouble;
/// Minimum value of t
double m_tMin;
/// Vector of detector data
std::vector<HistData> m_histData;
API::MatrixWorkspace_sptr squash(const API::MatrixWorkspace_sptr &ws,
const API::ITableWorkspace_sptr &phase,
const std::vector<double> &n0);
};
} // namespace Algorithms
......
Code/Mantid/Framework/Algorithms/src/PhaseQuadMuon.cpp
+ 172
− 389
View file @ 54de80c6
This diff is collapsed.
Code/Mantid/Framework/Algorithms/test/PhaseQuadMuonTest.h
+ 52
− 176
View file @ 54de80c6
......@@ -2,14 +2,10 @@
#define MANTID_ALGORITHMS_PHASEQUADMUONTEST_H_
#include <cxxtest/TestSuite.h>
#include "MantidDataHandling/LoadMuonNexus2.h"
#include "MantidAlgorithms/PhaseQuadMuon.h"
#include "MantidAPI/AlgorithmManager.h"
#include "MantidDataObjects/TableWorkspace.h"
#include "MantidAPI/TableRow.h"
#include <Poco/File.h>
#include <stdexcept>
using namespace Mantid::Algorithms;
using namespace Mantid::DataObjects;
using namespace Mantid::API;
......@@ -17,140 +13,43 @@ class PhaseQuadMuonTest : public CxxTest::TestSuite
{
public:
void testName()
void testTheBasics()
{
TS_ASSERT_EQUALS( phaseQuadMuon.name(), "PhaseQuad" );
IAlgorithm_sptr phaseQuad = AlgorithmManager::Instance().create("PhaseQuad");
TS_ASSERT_EQUALS(phaseQuad->name(), "PhaseQuad");
TS_ASSERT_EQUALS(phaseQuad->category(), "Muon");
TS_ASSERT_THROWS_NOTHING(phaseQuad->initialize());
TS_ASSERT(phaseQuad->isInitialized());
}
void testCategory()
{
TS_ASSERT_EQUALS( phaseQuadMuon.category(), "Muon" )
}
void testInit()
{
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.initialize() );
TS_ASSERT( phaseQuadMuon.isInitialized() )
}
void testExecPhaseList()
{
loader.initialize();
loader.setPropertyValue("Filename", "emu00006473.nxs");
loader.setPropertyValue("OutputWorkspace", "EMU6473");
TS_ASSERT_THROWS_NOTHING( loader.execute() );
TS_ASSERT_EQUALS(loader.isExecuted(),true);
MatrixWorkspace_sptr inputWs = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("EMU6473");
std::string filename("TestPhaseList.txt");
generatePhaseList(filename);
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("PhaseList", "TestPhaseList.txt") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("InputWorkspace", "EMU6473") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("OutputWorkspace", "EMU6473_out") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.execute() );
TS_ASSERT( phaseQuadMuon.isExecuted() );
MatrixWorkspace_sptr outputWs = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("EMU6473_out");
TS_ASSERT_EQUALS( outputWs->getNumberHistograms(), 2);
TS_ASSERT_EQUALS( outputWs->getSpectrum(0)->readX(), inputWs->getSpectrum(0)->readX() ); // Check outputWs X values
TS_ASSERT_EQUALS( outputWs->getSpectrum(1)->readX(), inputWs->getSpectrum(1)->readX() );
auto specReY = outputWs->getSpectrum(0)->readY();
auto specReE = outputWs->getSpectrum(0)->readE();
auto specImY = outputWs->getSpectrum(1)->readY();
auto specImE = outputWs->getSpectrum(1)->readE();
// Check real spectrum Y values
TS_ASSERT_DELTA ( specReY[ 0], -0.998265, 0.000001 );
TS_ASSERT_DELTA ( specReY[10], -0.997286, 0.000001 );
TS_ASSERT_DELTA ( specReY[20], -0.026196, 0.000001 );
TS_ASSERT_DELTA ( specReY[30], 0.017798, 0.000001 );
TS_ASSERT_DELTA ( specReY[40], 0.033196, 0.000001 );
TS_ASSERT_DELTA ( specReY[50], 0.025337, 0.000001 );
// Check real spectrum E values
TS_ASSERT_DELTA ( specReE[ 0], 135268, 1 );
TS_ASSERT_DELTA ( specReE[10], 145487, 1 );
TS_ASSERT_DELTA ( specReE[20], 0.00213851, 0.000001 );
TS_ASSERT_DELTA ( specReE[30], 0.00226644, 0.000001 );
TS_ASSERT_DELTA ( specReE[40], 0.00237071, 0.000001 );
TS_ASSERT_DELTA ( specReE[50], 0.00244977, 0.000001 );
// Check imaginary spectrum Y values
TS_ASSERT_DELTA ( specImY[ 0], -0.997455, 0.000001 );
TS_ASSERT_DELTA ( specImY[10], -0.993110, 0.000001 );
TS_ASSERT_DELTA ( specImY[20], 0.0099704, 0.000001 );
TS_ASSERT_DELTA ( specImY[30], 0.0300842, 0.000001 );
TS_ASSERT_DELTA ( specImY[40], 0.0285628, 0.000001 );
TS_ASSERT_DELTA ( specImY[50], 0.0300885, 0.000001 );
// Check imaginary spectrum E values
TS_ASSERT_DELTA ( specImE[ 0], 280312, 1 );
TS_ASSERT_DELTA ( specImE[10], 301487, 1 );
TS_ASSERT_DELTA ( specImE[20], 0.00316581, 0.000001 );
TS_ASSERT_DELTA ( specImE[30], 0.00332145, 0.000001 );
TS_ASSERT_DELTA ( specImE[40], 0.00343792, 0.000001 );
TS_ASSERT_DELTA ( specImE[50], 0.00357113, 0.000001 );
AnalysisDataService::Instance().remove("EMU6473"); // remove inputWs
AnalysisDataService::Instance().remove("EMU6473_out"); // remove OutputWs
Poco::File("TestPhaseList.txt").remove(); // remove phase list
}
void generatePhaseList (std::string filename)
{
std::ofstream ofile;
ofile.open(filename.c_str());
if (ofile.is_open())
{
// Write header
ofile << "MuSR 64 det 12705-12715" << std::endl;
ofile << "Top row of numbers are:" << std::endl;
ofile << "#histos, typ. first good bin#, typ. bin# when pulse over, mean lag." << std::endl;
ofile << "Tabulated numbers are, per histogram:" << std::endl;
ofile << "det ok, asymmetry, phase, lag, deadtime_c, deadtime_m." << std::endl;
ofile << "32 2 0 0" << std::endl;
// Write data
for (int i=0; i<16; i++)
{
ofile << "1 50.0 0.00 0.0 0.0 1" << std::endl;
ofile << "1 50.0 1.57 0.0 0.0 1" << std::endl;
}
ofile.close();
}
else
{
throw std::runtime_error("Unable to open file to write.");
}
}
void testExecPhaseTable()
{
loader.initialize();
loader.setPropertyValue("Filename", "emu00006473.nxs");
loader.setPropertyValue("OutputWorkspace", "EMU6473");
TS_ASSERT_THROWS_NOTHING( loader.execute() );
TS_ASSERT_EQUALS(loader.isExecuted(),true);
MatrixWorkspace_sptr inputWs = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("EMU6473");
// Load a muon dataset
IAlgorithm_sptr loader = AlgorithmManager::Instance().create("Load");
loader->setChild(true);
loader->initialize();
loader->setProperty("Filename", "emu00006473.nxs");
loader->setPropertyValue("OutputWorkspace", "out_ws");
loader->execute();
Workspace_sptr temp = loader->getProperty("OutputWorkspace");
MatrixWorkspace_sptr inputWs = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);
// Create and populate a detector table
boost::shared_ptr<ITableWorkspace>phaseTable(new Mantid::DataObjects::TableWorkspace);
generatePhaseTable(phaseTable);
AnalysisDataService::Instance().add("PhaseTable",phaseTable);
populatePhaseTable(phaseTable);
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("InputWorkspace", "EMU6473") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("OutputWorkspace", "EMU6473_out") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("PhaseTable", "PhaseTable") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("PulseOver", "60") );
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.setProperty("MeanLag", "0") );
// Run PhaseQuad
IAlgorithm_sptr phaseQuad = AlgorithmManager::Instance().create("PhaseQuad");
phaseQuad->setChild(true);
phaseQuad->initialize();
phaseQuad->setProperty("InputWorkspace",inputWs);
phaseQuad->setProperty("PhaseTable",phaseTable);
phaseQuad->setPropertyValue("OutputWorkspace", "out_ws");
TS_ASSERT_THROWS_NOTHING(phaseQuad->execute());
TS_ASSERT(phaseQuad->isExecuted());
TS_ASSERT_THROWS_NOTHING( phaseQuadMuon.execute() );
TS_ASSERT( phaseQuadMuon.isExecuted() );
MatrixWorkspace_sptr outputWs = AnalysisDataService::Instance().retrieveWS<MatrixWorkspace>("EMU6473_out");
// Get the output ws
MatrixWorkspace_sptr outputWs = phaseQuad->getProperty("OutputWorkspace");
TS_ASSERT_EQUALS( outputWs->getNumberHistograms(), 2);
TS_ASSERT_EQUALS( outputWs->getSpectrum(0)->readX(), inputWs->getSpectrum(0)->readX() ); // Check outputWs X values
......@@ -160,59 +59,36 @@ public:
auto specReE = outputWs->getSpectrum(0)->readE();
auto specImY = outputWs->getSpectrum(1)->readY();
auto specImE = outputWs->getSpectrum(1)->readE();
// Check real spectrum Y values
TS_ASSERT_DELTA ( specReY[ 0], -0.998265, 0.000001 );
TS_ASSERT_DELTA ( specReY[10], -0.997286, 0.000001 );
TS_ASSERT_DELTA ( specReY[20], -0.026196, 0.000001 );
TS_ASSERT_DELTA ( specReY[30], 0.017798, 0.000001 );
TS_ASSERT_DELTA ( specReY[40], 0.033196, 0.000001 );
TS_ASSERT_DELTA ( specReY[50], 0.025337, 0.000001 );
// Check real spectrum E values
TS_ASSERT_DELTA ( specReE[ 0], 135268, 1 );
TS_ASSERT_DELTA ( specReE[10], 145487, 1 );
TS_ASSERT_DELTA ( specReE[20], 0.00213851, 0.000001 );
TS_ASSERT_DELTA ( specReE[30], 0.00226644, 0.000001 );
TS_ASSERT_DELTA ( specReE[40], 0.00237071, 0.000001 );
TS_ASSERT_DELTA ( specReE[50], 0.00244977, 0.000001 );
// Check imaginary spectrum Y values
TS_ASSERT_DELTA ( specImY[ 0], -0.997455, 0.000001 );
TS_ASSERT_DELTA ( specImY[10], -0.993110, 0.000001 );
TS_ASSERT_DELTA ( specImY[20], 0.0099704, 0.000001 );
TS_ASSERT_DELTA ( specImY[30], 0.0300842, 0.000001 );
TS_ASSERT_DELTA ( specImY[40], 0.0285628, 0.000001 );
TS_ASSERT_DELTA ( specImY[50], 0.0300885, 0.000001 );
// Check imaginary spectrum E values
TS_ASSERT_DELTA ( specImE[ 0], 280312, 1 );
TS_ASSERT_DELTA ( specImE[10], 301487, 1 );
TS_ASSERT_DELTA ( specImE[20], 0.00316581, 0.000001 );
TS_ASSERT_DELTA ( specImE[30], 0.00332145, 0.000001 );
TS_ASSERT_DELTA ( specImE[40], 0.00343792, 0.000001 );
TS_ASSERT_DELTA ( specImE[50], 0.00357113, 0.000001 );
AnalysisDataService::Instance().remove("EMU6473"); // remove inputWs
AnalysisDataService::Instance().remove("EMU6473_out"); // remove OutputWs
AnalysisDataService::Instance().remove("PhaseTable"); // remove PhaseTable
// Check real Y values
TS_ASSERT_DELTA ( specReY[ 0], -0.9982, 0.0001 );
TS_ASSERT_DELTA ( specReY[20], -0.0251, 0.0001 );
TS_ASSERT_DELTA ( specReY[50], 0.0264, 0.0001 );
// Check real E values
TS_ASSERT_DELTA ( specReE[ 0], 0.0010, 0.0001);
TS_ASSERT_DELTA ( specReE[20], 0.0021, 0.0001 );
TS_ASSERT_DELTA ( specReE[50], 0.0024, 0.0001 );
// Check imaginary Y values
TS_ASSERT_DELTA ( specImY[ 0], -0.9974, 0.0001 );
TS_ASSERT_DELTA ( specImY[20], 0.0115, 0.0001 );
TS_ASSERT_DELTA ( specImY[50], 0.0316, 0.0001 );
// Check imaginary E values
TS_ASSERT_DELTA ( specImE[ 0], 0.0029, 0.0001 );
TS_ASSERT_DELTA ( specImE[20], 0.0031, 0.0001 );
TS_ASSERT_DELTA ( specImE[50], 0.0035, 0.0001 );
}
void generatePhaseTable (ITableWorkspace_sptr phaseTable)
void populatePhaseTable (ITableWorkspace_sptr phaseTable)
{
phaseTable->addColumn("bool","DetectorOK");
phaseTable->addColumn("double","DetectorAlpha");
phaseTable->addColumn("int","DetectorID");
phaseTable->addColumn("double", "DetectorAsymmetry");
phaseTable->addColumn("double","DetectorPhase");
phaseTable->addColumn("double","DetectorDeadTime");
for (int i=0; i<16; i++)
{
for (int i = 0; i < 16; i++) {
TableRow phaseRow1 = phaseTable->appendRow();
phaseRow1 << true << 50.0 << 0.0 << 0.0 ;
phaseRow1 << i << 1. << 0.;
TableRow phaseRow2 = phaseTable->appendRow();
phaseRow2 << true << 50.0 << 1.57 << 0.0 ;
phaseRow2 << i << 1. << 1.57;
}
}
private:
PhaseQuadMuon phaseQuadMuon;
Mantid::DataHandling::LoadMuonNexus2 loader;
};
#endif /* MANTID_ALGORITHMS_PHASEQUADMUONTEST_H_ */
\ No newline at end of file
Code/Mantid/docs/source/algorithms/PhaseQuad-v1.rst
+ 29
− 81
View file @ 54de80c6
......@@ -5,14 +5,14 @@
Overview
--------
This algorithm provides a method for combining signals from multiple detector segments to form two output
signals in quadrature phase. It is of particular use when working with muon spin rotation signals measured
using highly segmented detector arrays that are typically found at pulsed muon sources (ISIS has instruments
with up to 600 detector elements). The method allows information from individual detectors to be combined in
a way that makes full use of the dataset.
This algorithm is frequently run as a precursor to making a Rotating Reference Frame transformation of
the dataset using the algorithm :ref:`algm-RRFMuon`. Both algorithms are fully described in the article
This algorithm provides a method for combining signals from multiple detector segments to form two output
signals in quadrature phase. It is of particular use when working with muon spin rotation signals measured
using highly segmented detector arrays that are typically found at pulsed muon sources (ISIS has instruments
with up to 600 detector elements). The method allows information from individual detectors to be combined in
a way that makes full use of the dataset.
This algorithm is frequently run as a precursor to making a Rotating Reference Frame transformation of
the dataset using the algorithm :ref:`algm-RRFMuon`. Both algorithms are fully described in the article
by T.M. Riseman and J.H. Brewer [Hyp. Int., 65, (1990), 1107].
......@@ -23,96 +23,44 @@ by T.M. Riseman and J.H. Brewer [Hyp. Int., 65, (1990), 1107].
Description
-----------
Assuming that the *InputWorkspace* contains measured
counts as a function of time, the algorithm returns a workspace
Assuming that the *InputWorkspace* contains measured counts as a function of time,
and *PhaseTable* contains the detector phases and asymmetries, the algorithm returns a workspace
containing two spectra (squashograms) as a function of the same time binning.
Either *PhaseList* or *PhaseTable* must be provided. Note: The first column is ignored in both cases.
*PhaseList* is expected to have a six-row header, including overall information, and 6 columns:
1. Boolean type, containing if detector is OK
2. Double type, containing detector asymmetry
3. Double type, containing detector phase
4. Double type, containing detector lag
5. Double type, containing detector deadtime :math:`c`
6. Double type, containing detector deadtime :math:`m`
*PhaseTable* is expected to have four columns:
1. Boolean type, containing if detector is OK
2. Double type, containing detector asymmetry
3. Double type, containing detector phase
4. Double type, containing detector deadtime
*PhaseTable* is expected to have three columns, corresponding to the detector id, its asymmetry
and its phase.
Usage
-----
.. include:: ../usagedata-note.txt
**Example - Computing squashograms from PhaseList:**
.. testcode:: ExPhaseQuadList
# Load a set of spectra from a EMU file
ws = LoadMuonNexus('emu00006473.nxs')
# Create a PhaseList file with some arbitrary detector information
import os
phaselist_path = os.path.join(os.path.expanduser("~"),"PhaseList.txt")
phaselist_file = open(phaselist_path,'w')
phaselist_file.write("MuSR\n")
phaselist_file.write("Dummy line\n")
phaselist_file.write("Dummy line\n")
phaselist_file.write("Dummy line\n")
phaselist_file.write("Dummy line\n")
phaselist_file.write("32 0 60 0.0\n")
for i in range(0,16):
phaselist_file.write("1 50.0 0.00 0 0 1\n")
phaselist_file.write("1 50.0 1.57 0 0 1\n")
phaselist_file.close()
ows = PhaseQuad(InputWorkspace='ws',PhaseList=phaselist_path)
print "Output workspace contains", ows.getNumberHistograms(), "histograms"
**Example - Computing squashograms:**
Output:
.. testoutput:: ExPhaseQuadList
Output workspace contains 2 histograms
.. testcleanup:: ExPhaseQuadList
.. testcode:: ExPhaseQuad
import os
try:
os.remove(phaselist_path)
except OSError:
pass
**Example - Computing squashograms from PhaseTable:**
.. testcode:: ExPhaseQuadTable
# Load a set of spectra from a EMU file
ws = LoadMuonNexus('emu00006473.nxs')
from math import pi
Load(Filename='MUSR00022725.nxs', OutputWorkspace='MUSR00022725')
CropWorkspace(InputWorkspace='MUSR00022725', OutputWorkspace='MUSR00022725', XMin=0, EndWorkspaceIndex=63)
# Create a PhaseTable with some arbitrary detector information
tab = CreateEmptyTableWorkspace()
tab.addColumn('bool', 'Status')
tab.addColumn('double', 'Asymmetry')
tab.addColumn('int', 'DetID')
tab.addColumn('double', 'Phase')
tab.addColumn('double', 'DeadTime')
for i in range(0,16):
tab.addRow([1, 50.0, 0.00, 0])
tab.addRow([1, 50.0, 1.57, 0])
ows = PhaseQuad(InputWorkspace='ws',PhaseTable='tab')
print "Output workspace contains", ows.getNumberHistograms(), "histograms"
tab.addColumn('double', 'Asym')
for i in range(0,32):
phi = 2*pi*i/32.
tab.addRow([1, 0.2, phi])
for i in range(0,32):
phi = 2*pi*i/32.
tab.addRow([1, 0.2, phi])
ows = PhaseQuad(InputWorkspace='MUSR00022725', PhaseTable='tab')
print "Output workspace has", ows.getNumberHistograms(), "histograms"
Output:
.. testoutput:: ExPhaseQuadTable
.. testoutput:: ExPhaseQuad
Output workspace contains 2 histograms
Output workspace has 2 histograms
.. categories::
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment