-
Notifications
You must be signed in to change notification settings - Fork 121
/
AlphaCalc.cpp
149 lines (120 loc) · 5.16 KB
/
AlphaCalc.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
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidAlgorithms/AlphaCalc.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidKernel/ArrayProperty.h"
#include <cmath>
#include <vector>
namespace Mantid {
namespace Algorithms {
using namespace Kernel;
using API::Progress;
// Register the class into the algorithm factory
DECLARE_ALGORITHM(AlphaCalc)
/** Initialisation method. Declares properties to be used in algorithm.
*
*/
void AlphaCalc::init() {
declareProperty(Kernel::make_unique<API::WorkspaceProperty<>>(
"InputWorkspace", "", Direction::Input),
"Name of the input workspace");
std::vector<int> forwardDefault{1};
declareProperty(
Kernel::make_unique<ArrayProperty<int>>("ForwardSpectra", forwardDefault),
"The spectra numbers of the forward group (default to 1)");
std::vector<int> backwardDefault{2};
declareProperty(Kernel::make_unique<ArrayProperty<int>>("BackwardSpectra",
backwardDefault),
"The spectra numbers of the backward group (default to 2)");
declareProperty("FirstGoodValue", EMPTY_DBL(),
"First good value (default lowest value of x)",
Direction::Input);
declareProperty("LastGoodValue", EMPTY_DBL(),
"Last good value (default highest value of x)",
Direction::Input);
declareProperty("Alpha", 1.0, "The alpha efficiency (default to 1.0)",
Direction::Output);
}
/** Executes the algorithm
*
*/
void AlphaCalc::exec() {
std::vector<int> forwardSpectraList = getProperty("ForwardSpectra");
std::vector<int> backwardSpectraList = getProperty("BackwardSpectra");
double alpha = getProperty("Alpha");
// no point in attempting to calculate alpha if input workspace on contains
// one spectra
API::MatrixWorkspace_sptr inputWS = getProperty("InputWorkspace");
if (inputWS->getNumberHistograms() < 2) {
g_log.error()
<< "Can't calculate alpha value for workspace which"
<< " contains one spectrum. A default value of alpha = 1.0 is returned";
setProperty("Alpha", alpha);
return;
}
// if for some reason the size of forward and backward lists are zero
// default these to their defaults
if (forwardSpectraList.empty())
forwardSpectraList.push_back(1);
if (backwardSpectraList.empty())
backwardSpectraList.push_back(2);
// first step is to create two workspaces which groups all forward and
// backward spectra
API::IAlgorithm_sptr groupForward = createChildAlgorithm("GroupDetectors");
groupForward->setProperty("InputWorkspace", inputWS);
groupForward->setProperty("OutputWorkspace", "tmp");
groupForward->setProperty("SpectraList", forwardSpectraList);
groupForward->setProperty("KeepUngroupedSpectra", false);
groupForward->execute();
API::MatrixWorkspace_sptr forwardWS =
groupForward->getProperty("OutputWorkspace");
API::IAlgorithm_sptr groupBackward = createChildAlgorithm("GroupDetectors");
groupBackward->setProperty("InputWorkspace", inputWS);
groupBackward->setProperty("OutputWorkspace", "tmp");
groupBackward->setProperty("SpectraList", backwardSpectraList);
groupBackward->setProperty("KeepUngroupedSpectra", false);
groupBackward->execute();
API::MatrixWorkspace_sptr backwardWS =
groupBackward->getProperty("OutputWorkspace");
// calculate sum of forward counts
double firstGoodvalue = getProperty("FirstGoodValue");
double lastGoodvalue = getProperty("LastGoodValue");
API::IAlgorithm_sptr integr = createChildAlgorithm("Integration");
integr->setProperty("InputWorkspace", forwardWS);
integr->setPropertyValue("OutputWorkspace", "tmp");
if (firstGoodvalue != EMPTY_DBL())
integr->setProperty("RangeLower", firstGoodvalue);
if (lastGoodvalue != EMPTY_DBL())
integr->setProperty("RangeUpper", lastGoodvalue);
integr->execute();
API::MatrixWorkspace_sptr out = integr->getProperty("OutputWorkspace");
double sumForward = out->readY(0)[0];
if (sumForward < 0) {
g_log.error() << "Sum of forward detector counts is negative."
<< "Therefore can't calculate alpha. Return alpha = 1.0.";
setProperty("Alpha", alpha);
return;
}
// calculate sum of backward counts
API::IAlgorithm_sptr integrB = createChildAlgorithm("Integration");
integrB->setProperty("InputWorkspace", backwardWS);
integrB->setPropertyValue("OutputWorkspace", "tmp");
if (firstGoodvalue != EMPTY_DBL())
integrB->setProperty("RangeLower", firstGoodvalue);
if (lastGoodvalue != EMPTY_DBL())
integrB->setProperty("RangeUpper", lastGoodvalue);
integrB->execute();
out = integrB->getProperty("OutputWorkspace");
double sumBackward = out->readY(0)[0];
if (sumBackward <= 0) {
g_log.error() << "Sum of backward detector counts is negative or zero."
<< "Therefore can't calculate alpha. Return alpha = 1.0.";
setProperty("Alpha", alpha);
return;
}
// finally calculate alpha
setProperty("Alpha", sumForward / sumBackward);
}
} // namespace Algorithm
} // namespace Mantid