-
Notifications
You must be signed in to change notification settings - Fork 122
/
NormaliseVanadium.cpp
133 lines (109 loc) · 4.43 KB
/
NormaliseVanadium.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
// Mantid Repository : https://github.com/mantidproject/mantid
//
// Copyright © 2018 ISIS Rutherford Appleton Laboratory UKRI,
// NScD Oak Ridge National Laboratory, European Spallation Source,
// Institut Laue - Langevin & CSNS, Institute of High Energy Physics, CAS
// SPDX - License - Identifier: GPL - 3.0 +
#include "MantidCrystal/NormaliseVanadium.h"
#include "MantidAPI/Axis.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/SpectrumInfo.h"
#include "MantidAPI/WorkspaceFactory.h"
#include "MantidGeometry/Instrument.h"
#include "MantidKernel/BoundedValidator.h"
#include "MantidKernel/Fast_Exponential.h"
#include "MantidKernel/Unit.h"
/* Following A.J.Schultz's anvred, scaling the vanadium spectra:
*/
namespace Mantid {
namespace Crystal {
// Register the class into the algorithm factory
DECLARE_ALGORITHM(NormaliseVanadium)
using namespace Kernel;
using namespace Geometry;
using namespace API;
using namespace DataObjects;
NormaliseVanadium::NormaliseVanadium() : API::Algorithm() {}
void NormaliseVanadium::init() {
// The input workspace must have an instrument and units of wavelength
auto wsValidator = std::make_shared<InstrumentValidator>();
declareProperty(std::make_unique<WorkspaceProperty<>>(
"InputWorkspace", "", Direction::Input, wsValidator),
"The X values for the input workspace must be in units of "
"wavelength or TOF");
declareProperty(std::make_unique<WorkspaceProperty<>>("OutputWorkspace", "",
Direction::Output),
"Output workspace name");
auto mustBePositive = std::make_shared<BoundedValidator<double>>();
mustBePositive->setLower(0.0);
declareProperty("Wavelength", 1.0, mustBePositive,
"Normalizes spectra to this wavelength");
}
void NormaliseVanadium::exec() {
// Retrieve the input workspace
m_inputWS = getProperty("InputWorkspace");
std::string unitStr = m_inputWS->getAxis(0)->unit()->unitID();
// Get the input parameters
double lambdanorm = getProperty("Wavelength"); // in 1/cm
MatrixWorkspace_sptr correctionFactors =
WorkspaceFactory::Instance().create(m_inputWS);
const auto numHists = static_cast<int64_t>(m_inputWS->getNumberHistograms());
const auto specSize = static_cast<int64_t>(m_inputWS->blocksize());
// If sample not at origin, shift cached positions.
const auto &spectrumInfo = m_inputWS->spectrumInfo();
double L1 = spectrumInfo.l1();
Progress prog(this, 0.0, 1.0, numHists);
// Loop over the spectra
PARALLEL_FOR_IF(Kernel::threadSafe(*m_inputWS, *correctionFactors))
for (int64_t i = 0; i < int64_t(numHists); ++i) {
// PARALLEL_START_INTERUPT_REGION //FIXME: Restore
// Get a reference to the Y's in the output WS for storing the factors
auto &Y = correctionFactors->mutableY(i);
auto &E = correctionFactors->mutableE(i);
// Copy over bin boundaries
const auto &inSpec = m_inputWS->getSpectrum(i);
correctionFactors->setSharedX(i, inSpec.sharedX());
const auto &Yin = inSpec.y();
const auto &Ein = inSpec.e();
// If no detector is found, skip onto the next spectrum
if (!spectrumInfo.hasDetectors(i))
continue;
// This is the scattered beam direction
double L2 = spectrumInfo.l2(i);
// Two-theta = polar angle = scattering angle = between +Z vector and the
// scattered beam
double scattering = spectrumInfo.twoTheta(i);
Mantid::Kernel::Units::Wavelength wl;
auto timeflight = inSpec.points();
if (unitStr == "TOF")
wl.fromTOF(timeflight.mutableRawData(), timeflight.mutableRawData(), L1,
L2, scattering, 0, 0, 0);
// Loop through the bins in the current spectrum
double lambp = 0;
double lambm = 0;
double normp = 0;
double normm = 0;
for (int64_t j = 0; j < specSize; j++) {
const double lambda = timeflight[j];
if (lambda > lambdanorm) {
lambp = lambda;
normp = Yin[j];
break;
}
lambm = lambda;
normm = Yin[j];
}
double normvalue =
normm + (lambdanorm - lambm) * (normp - normm) / (lambp - lambm);
for (int64_t j = 0; j < specSize; j++) {
Y[j] = Yin[j] / normvalue;
E[j] = Ein[j] / normvalue;
}
prog.report();
// PARALLEL_END_INTERUPT_REGION
}
// PARALLEL_CHECK_INTERUPT_REGION
setProperty("OutputWorkspace", correctionFactors);
}
} // namespace Crystal
} // namespace Mantid