-
Notifications
You must be signed in to change notification settings - Fork 122
/
ApplyDeadTimeCorr.cpp
152 lines (127 loc) · 5.44 KB
/
ApplyDeadTimeCorr.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
/*WIKI*
Apply deadtime correction to each spectra of a workspace. Define:
<math>{\displaystyle{N}}</math> = true count
<math>{\displaystyle{M}}</math> = measured count
<math>{\displaystyle{t_{dead}}}</math> = dead-time
<math>{\displaystyle{t_{bin}}}</math> = time bin width
<math>{\displaystyle{F}}</math> = Number of good frames
Then this algorithm assumes that the InputWorkspace contains measured counts as a
function of TOF and returns a workspace containing true counts as a function of the
same TOF binning according to
:<math> N = \frac{M}{(1-M*(\frac{t_{dead}}{t_{bin}*F}))} </math>
*WIKI*/
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/ApplyDeadTimeCorr.h"
#include "MantidKernel/System.h"
#include "MantidKernel/TimeSeriesProperty.h"
#include "MantidKernel/PropertyWithValue.h"
#include "boost/lexical_cast.hpp"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidAPI/ITableWorkspace.h"
#include "MantidAPI/TableRow.h"
#include <iostream>
#include <cmath>
using std::string;
using namespace Mantid::Kernel;
using namespace Mantid::API;
namespace Mantid
{
namespace Algorithms
{
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(ApplyDeadTimeCorr)
//----------------------------------------------------------------------------------------------
/** Initialize the algorithm's properties.
*/
void ApplyDeadTimeCorr::init()
{
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::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");
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") )
{
double numGoodFrames = boost::lexical_cast<double>(run.getProperty("goodfrm")->value());
if (numGoodFrames == 0)
{
throw std::runtime_error("Number of good frames in the workspace is zero");
}
// 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
{
// 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
{
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() << "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