-
Notifications
You must be signed in to change notification settings - Fork 122
/
EQSANSPatchSensitivity.cpp
137 lines (118 loc) · 4.73 KB
/
EQSANSPatchSensitivity.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
#include "MantidWorkflowAlgorithms/EQSANSPatchSensitivity.h"
#include "MantidAPI/MatrixWorkspace.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/cow_ptr.h"
namespace Mantid {
namespace WorkflowAlgorithms {
// Register the algorithm into the AlgorithmFactory
DECLARE_ALGORITHM(EQSANSPatchSensitivity)
using namespace Kernel;
using namespace API;
using namespace Geometry;
void EQSANSPatchSensitivity::init() {
declareProperty(
make_unique<WorkspaceProperty<>>("Workspace", "", Direction::InOut),
"Input sensitivity workspace to be patched");
declareProperty(
make_unique<WorkspaceProperty<>>("PatchWorkspace", "", Direction::Input),
"Workspace defining the patch. Masked detectors will be patched.");
declareProperty("UseLinearRegression", true, "If true, a linear regression "
"will be used instead of "
"computing the average");
declareProperty("OutputMessage", "", Direction::Output);
}
void EQSANSPatchSensitivity::exec() {
MatrixWorkspace_sptr inputWS = getProperty("Workspace");
MatrixWorkspace_sptr patchWS = getProperty("PatchWorkspace");
bool useRegression = getProperty("UseLinearRegression");
const int nx_pixels = static_cast<int>(
inputWS->getInstrument()->getNumberParameter("number-of-x-pixels")[0]);
const int ny_pixels = static_cast<int>(
inputWS->getInstrument()->getNumberParameter("number-of-y-pixels")[0]);
const int numberOfSpectra = static_cast<int>(inputWS->getNumberHistograms());
const auto &spectrumInfo = patchWS->spectrumInfo();
const auto &inSpectrumInfo = inputWS->spectrumInfo();
// Loop over all tubes and patch as necessary
for (int i = 0; i < nx_pixels; i++) {
std::vector<int> patched_ids;
int nUnmasked = 0;
double totalUnmasked = 0.0;
double errorUnmasked = 0.0;
double sumXY = 0.0;
double sumX = 0.0;
double sumX2 = 0.0;
double sumY = 0.0;
progress(0.9 * i / nx_pixels, "Processing patch");
for (int j = 0; j < ny_pixels; j++) {
// EQSANS-specific: get detector ID from pixel coordinates
int iDet = ny_pixels * i + j;
if (iDet > numberOfSpectra) {
g_log.notice() << "Got an invalid detector ID " << iDet << '\n';
continue;
}
// If this detector is a monitor, skip to the next one
if (spectrumInfo.isMonitor(iDet))
continue;
const MantidVec &YValues = inputWS->readY(iDet);
const MantidVec &YErrors = inputWS->readE(iDet);
// If this detector is masked, skip to the next one
if (spectrumInfo.isMasked(i))
patched_ids.push_back(iDet);
else {
if (!inSpectrumInfo.isMasked(iDet)) {
double yPosition = spectrumInfo.position(iDet).Y();
totalUnmasked += YErrors[0] * YErrors[0] * YValues[0];
errorUnmasked += YErrors[0] * YErrors[0];
nUnmasked++;
sumXY += yPosition * YValues[0];
sumX += yPosition;
sumX2 += yPosition * yPosition;
sumY += YValues[0];
}
}
}
if (nUnmasked > 0 && errorUnmasked > 0) {
sumXY /= nUnmasked;
sumX /= nUnmasked;
sumX2 /= nUnmasked;
sumY /= nUnmasked;
double beta = (sumXY - sumX * sumY) / (sumX2 - sumX * sumX);
double alpha = sumY - beta * sumX;
double error = sqrt(errorUnmasked) / nUnmasked;
double average = totalUnmasked / errorUnmasked;
// Apply patch
progress(0.91, "Applying patch");
auto &spectrumInfo = inputWS->mutableSpectrumInfo();
for (auto patched_id : patched_ids) {
if (!spectrumInfo.hasDetectors(patched_id)) {
g_log.warning() << "Spectrum " << patched_id
<< " has no detector, skipping (not clearing mask)\n";
continue;
}
MantidVec &YValues = inputWS->dataY(patched_id);
MantidVec &YErrors = inputWS->dataE(patched_id);
if (useRegression) {
YValues[0] = alpha + beta * spectrumInfo.position(patched_id).Y();
YErrors[0] = error;
} else {
YValues[0] = average;
YErrors[0] = error;
}
spectrumInfo.setMasked(patched_id, false);
}
}
}
// Call Calculate efficiency to renormalize
progress(0.91, "Renormalizing");
IAlgorithm_sptr effAlg = createChildAlgorithm("CalculateEfficiency");
effAlg->setProperty("InputWorkspace", inputWS);
effAlg->setProperty("OutputWorkspace", inputWS);
effAlg->execute();
inputWS = effAlg->getProperty("OutputWorkspace");
setProperty("Workspace", inputWS);
setProperty("OutputMessage",
"Applied wavelength-dependent sensitivity correction");
}
} // namespace WorkflowAlgorithms
} // namespace Mantid