Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/feature/7229_muon_plotbylog_dead…
Browse files Browse the repository at this point in the history 
…time'
  • Loading branch information
@Anders-Markvardsen
Anders-Markvardsen committed Oct 16, 2013
2 parents 7b82e2c + a58edd5 commit a77a209
Show file tree
Hide file tree
Showing 11 changed files with 1,206 additions and 333 deletions.
36 changes: 30 additions & 6 deletions Code/Mantid/Framework/Algorithms/inc/MantidAlgorithms/PlotAsymmetryByLogValue.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
// Includes
//----------------------------------------------------------------------
#include "MantidAPI/Algorithm.h"
#include "MantidDataObjects/Workspace2D.h"

namespace Mantid
{
Expand All @@ -18,6 +19,8 @@ namespace Mantid

namespace Algorithms
{
using namespace API;
using namespace DataObjects;
/**Takes a muon workspace as input and sums all the spectra into two spectra which represent
the two detector groupings. The resultant spectra are used to calculate (F-aB) / (F+aB) the results of which
are stored in the output workspace.
Expand Down Expand Up @@ -59,7 +62,7 @@ namespace Mantid
{
public:
/// Default constructor
PlotAsymmetryByLogValue() : API::Algorithm() {};
PlotAsymmetryByLogValue() : Algorithm() {};
/// Destructor
virtual ~PlotAsymmetryByLogValue() {};
/// Algorithm's name for identification overriding a virtual method
Expand All @@ -77,15 +80,36 @@ namespace Mantid
void exec();

/// Calculate the integral asymmetry for a workspace (single period)
void calcIntAsymmetry(API::MatrixWorkspace_sptr ws, double& Y, double& E);
void calcIntAsymmetry(MatrixWorkspace_sptr ws, double& Y, double& E);

/// Calculate the integral asymmetry for a workspace (red & green)
void calcIntAsymmetry(API::MatrixWorkspace_sptr ws_red,
API::MatrixWorkspace_sptr ws_geen,double& Y, double& E);
void calcIntAsymmetry(MatrixWorkspace_sptr ws_red, MatrixWorkspace_sptr ws_geen,double& Y, double& E);
/// Group detectors
void groupDetectors(API::MatrixWorkspace_sptr& ws,const std::vector<int>& spectraList);
void groupDetectors(MatrixWorkspace_sptr& ws,const std::vector<int>& spectraList);
/// Get log value
double getLogValue(API::MatrixWorkspace& ws,const std::string& logName);
double getLogValue(MatrixWorkspace& ws,const std::string& logName);

/// Create a Dead Time Table using given list of dead times
Workspace_sptr deadTimesToTable(const std::vector<double>& deadTimes, size_t numDetectors = 0);

/// Creates Dead Time Table using all the data between begin and end
ITableWorkspace_sptr createDeadTimeTable( std::vector<double>::const_iterator begin,
std::vector<double>::const_iterator end);

/// Runs an appropriate applyDeadTimeCorrection function depending on the type of workspaces
Workspace_sptr applyDeadTimeCorrection(Workspace_sptr deadTimeWs, Workspace_sptr ws);

/// Applies DTC to a group of workspaces using a single table
WorkspaceGroup_sptr applyDeadTimeCorrection(ITableWorkspace_sptr deadTimeTable,
WorkspaceGroup_sptr wsGroup);

/// Applies DTC to a group of workspace using a group of tables
WorkspaceGroup_sptr applyDeadTimeCorrection(WorkspaceGroup_sptr deadTimeGroup,
WorkspaceGroup_sptr wsGroup);

/// Runs ApplyDeadTimeCorr to apply DTC to a workspace using a table
Workspace2D_sptr applyDeadTimeCorrection(ITableWorkspace_sptr deadTimeTable,
Workspace2D_sptr ws);

/// Stores property "Int"
bool m_int;
Expand Down
144 changes: 71 additions & 73 deletions Code/Mantid/Framework/Algorithms/src/ApplyDeadTimeCorr.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -46,104 +46,102 @@ DECLARE_ALGORITHM(ApplyDeadTimeCorr)
*/
void ApplyDeadTimeCorr::init()
{
this->setWikiSummary("Apply deadtime correction to each spectra of a workspace.");
this->setWikiSummary("Apply deadtime correction to each spectra of a workspace.");

declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "",
Direction::Input), "The name of the input workspace containing measured counts");
declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("InputWorkspace", "",
Direction::Input), "The name of the input workspace containing measured counts");

declareProperty(new API::WorkspaceProperty<API::ITableWorkspace>("DeadTimeTable", "",
Direction::Input), "Name of the Dead Time Table");
declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace","",Direction::Output),
"The name of the output workspace containing corrected counts");
declareProperty(new API::WorkspaceProperty<API::ITableWorkspace>("DeadTimeTable", "",
Direction::Input), "Name of the Dead Time Table");

declareProperty(new API::WorkspaceProperty<API::MatrixWorkspace>("OutputWorkspace","",
Direction::Output), "The name of the output workspace containing corrected counts");
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
void ApplyDeadTimeCorr::exec()
{
// Get pointers to the workspace and dead time table
API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
API::ITableWorkspace_sptr deadTimeTable = getProperty("DeadTimeTable");
// Get pointers to the workspace and dead time table
API::MatrixWorkspace_sptr inputWs = getProperty("InputWorkspace");
API::ITableWorkspace_sptr deadTimeTable = getProperty("DeadTimeTable");

if (!(deadTimeTable->rowCount() > inputWs->getNumberHistograms() ) )
{
// Get number of good frames from Run object. This also serves as
// a test to see if valid input workspace has been provided
if (!(deadTimeTable->rowCount() > inputWs->getNumberHistograms() ) )
{
// Get number of good frames from Run object. This also serves as
// a test to see if valid input workspace has been provided

const API::Run & run = inputWs->run();
if ( run.hasProperty("goodfrm") )
const API::Run & run = inputWs->run();
if ( run.hasProperty("goodfrm") )
{
double numGoodFrames = boost::lexical_cast<double>(run.getProperty("goodfrm")->value());

// Duplicate the input workspace. Only need to change Y values based on dead time corrections
IAlgorithm_sptr duplicate = createChildAlgorithm("CloneWorkspace");
duplicate->initialize();
duplicate->setProperty<Workspace_sptr>("InputWorkspace", boost::dynamic_pointer_cast<Workspace>(inputWs));
duplicate->execute();
Workspace_sptr temp = duplicate->getProperty("OutputWorkspace");
MatrixWorkspace_sptr outputWs = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);

// Presumed to be the same for all data
double timeBinWidth(inputWs->dataX(0)[1] - inputWs->dataX(0)[0]);

if (timeBinWidth != 0)
{
try
{
double numGoodFrames = boost::lexical_cast<double>(run.getProperty("goodfrm")->value());

// Duplicate the input workspace. Only need to change Y values based on dead time corrections
IAlgorithm_sptr duplicate = createChildAlgorithm("CloneWorkspace");
duplicate->initialize();
duplicate->setProperty<Workspace_sptr>("InputWorkspace", boost::dynamic_pointer_cast<Workspace>(inputWs));
duplicate->execute();
Workspace_sptr temp = duplicate->getProperty("OutputWorkspace");
MatrixWorkspace_sptr outputWs = boost::dynamic_pointer_cast<MatrixWorkspace>(temp);
// Apply Dead Time
for (size_t i=0; i<deadTimeTable->rowCount(); ++i)
{
API::TableRow deadTimeRow = deadTimeTable->getRow(i);
size_t index = static_cast<size_t>(inputWs->getIndexFromSpectrumNumber(static_cast<int>(deadTimeRow.Int(0) ) ) );

// Presumed to be the same for all data
double timeBinWidth(inputWs->dataX(0)[1] - inputWs->dataX(0)[0]);

if (timeBinWidth != 0)
{
try
{
// Apply Dead Time
for (size_t i=0; i<deadTimeTable->rowCount(); ++i)
{
API::TableRow deadTimeRow = deadTimeTable->getRow(i);
size_t index = static_cast<size_t>(inputWs->getIndexFromSpectrumNumber(static_cast<int>(deadTimeRow.Int(0) ) ) );

for(size_t j=0; j<inputWs->blocksize(); ++j)
{
double temp(1 - inputWs->dataY(index)[j] * (deadTimeRow.Double(1)/ (timeBinWidth * numGoodFrames) ) );
if (temp != 0)
{
outputWs->dataY(index)[j] = inputWs->dataY(index)[j] / temp;
}
else
{
g_log.error() << "1 - MeasuredCount * (Deadtime/TimeBin width is currently (" << temp << "). Can't divide by this amount." << std::endl;

throw std::invalid_argument("Can't divide by 0");
}
}
}

setProperty("OutputWorkspace", outputWs);
}
catch(std::runtime_error&)
{
throw std::invalid_argument("Invalid argument for algorithm.");
}
}
else
for(size_t j=0; j<inputWs->blocksize(); ++j)
{
g_log.error() << "The time bin width is currently (" << timeBinWidth << "). Can't divide by this amount." << std::endl;
double temp(1 - inputWs->dataY(index)[j] * (deadTimeRow.Double(1)/ (timeBinWidth * numGoodFrames) ) );
if (temp != 0)
{
outputWs->dataY(index)[j] = inputWs->dataY(index)[j] / temp;
}
else
{
g_log.error() << "1 - MeasuredCount * (Deadtime/TimeBin width is currently (" << temp << "). Can't divide by this amount." << std::endl;

throw std::invalid_argument("Can't divide by 0");
}
}
}

setProperty("OutputWorkspace", outputWs);
}
else
catch(std::runtime_error&)
{
g_log.error() << "To calculate Muon deadtime requires that goodfrm (number of good frames) "
<< "is stored in InputWorkspace Run object\n";
throw std::invalid_argument("Invalid argument for algorithm.");
}
}
else
{
g_log.error() << "The time bin width is currently (" << timeBinWidth << "). Can't divide by this amount." << std::endl;

throw std::invalid_argument("Can't divide by 0");
}
}
else
{
g_log.error() << "Row count(" << deadTimeTable->rowCount() << ") of Dead time table is bigger than the Number of Histograms("
<< inputWs->getNumberHistograms() << ")." << std::endl;

throw std::invalid_argument("Row count was bigger than the Number of Histograms.");
g_log.error() << "To calculate Muon deadtime requires that goodfrm (number of good frames) "
<< "is stored in InputWorkspace Run object\n";
}
}
else
{
g_log.error() << "Row count(" << deadTimeTable->rowCount() << ") of Dead time table is bigger than the Number of Histograms("
<< inputWs->getNumberHistograms() << ")." << std::endl;

throw std::invalid_argument("Row count was bigger than the Number of Histograms.");
}
}



} // namespace Mantid
} // namespace Algorithms

0 comments on commit a77a209

Please sign in to comment.