DiffRotDiscreteCircle.cpp 8.46 KB
Newer Older
1
2
3
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
4
#include "MantidCurveFitting/Functions/DiffRotDiscreteCircle.h"
5

6
#include "MantidAPI/FunctionFactory.h"
7
#include "MantidAPI/MatrixWorkspace.h"
8
#include "MantidAPI/ParameterTie.h"
9
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
10
#include "MantidGeometry/IDetector.h"
11
#include "MantidKernel/Exception.h"
12
13
14
15
16
17
18
19
#include "MantidKernel/UnitConversion.h"

#include <cmath>
#include <limits>

namespace {
Mantid::Kernel::Logger g_log("DiffSphere");
}
20

21
22
namespace Mantid {
namespace CurveFitting {
23
namespace Functions {
Nick Draper's avatar
Nick Draper committed
24

25
using namespace CurveFitting;
26
using namespace Constraints;
27

28
using namespace API;
29

30
31
using namespace Geometry;

32
33
34
DECLARE_FUNCTION(ElasticDiffRotDiscreteCircle)
DECLARE_FUNCTION(InelasticDiffRotDiscreteCircle)
DECLARE_FUNCTION(DiffRotDiscreteCircle)
35

36
37
38
39
40
41
ElasticDiffRotDiscreteCircle::ElasticDiffRotDiscreteCircle() {
  // declareParameter("Height", 1.0); //parameter "Height" already declared in
  // constructor of base class DeltaFunction
  declareParameter("Radius", 1.0, "Circle radius [Angstroms] ");
  declareAttribute("Q", API::IFunction::Attribute(0.5));
  declareAttribute("N", API::IFunction::Attribute(3));
42
}
43

44
void ElasticDiffRotDiscreteCircle::init() {
45
  // Ensure positive values for Height and Radius
46
47
48
  BoundaryConstraint *HeightConstraint = new BoundaryConstraint(
      this, "Height", std::numeric_limits<double>::epsilon(), true);
  addConstraint(HeightConstraint);
49

50
51
52
  BoundaryConstraint *RadiusConstraint = new BoundaryConstraint(
      this, "Radius", std::numeric_limits<double>::epsilon(), true);
  addConstraint(RadiusConstraint);
53
54
}

55
56
57
58
double ElasticDiffRotDiscreteCircle::HeightPrefactor() const {
  const double R = getParameter("Radius");
  const double Q = getAttribute("Q").asDouble();
  const int N = getAttribute("N").asInt();
59
  double aN = 0;
60
61
62
  for (int k = 1; k < N; k++) {
    double x = 2 * Q * R * sin(M_PI * k / N);
    aN += sin(x) / x; // spherical Besell function of order zero j0==sin(x)/x
63
64
  }
  aN += 1; // limit for j0 when k==N, or x==0
65
  return aN / N;
66
67
}

68
69
70
71
72
73
74
InelasticDiffRotDiscreteCircle::InelasticDiffRotDiscreteCircle()
    : m_hbar(0.658211626) {
  declareParameter("Intensity", 1.0, "scaling factor [arbitrary units]");
  declareParameter("Radius", 1.0, "Circle radius [Angstroms]");
  declareParameter("Decay", 1.0, "Inverse of transition rate, in nanoseconds "
                                 "if energy in micro-ev, or picoseconds if "
                                 "energy in mili-eV");
75
  declareParameter("Shift", 0.0, "Shift in domain");
76

Dan Nixon's avatar
Dan Nixon committed
77
  declareAttribute("Q", API::IFunction::Attribute(EMPTY_DBL()));
78
  declareAttribute("WorkspaceIndex", API::IFunction::Attribute(0));
79
  declareAttribute("N", API::IFunction::Attribute(3));
80
}
81

82
void InelasticDiffRotDiscreteCircle::init() {
83
  // Ensure positive values for Intensity, Radius, and decay
84
85
  BoundaryConstraint *IntensityConstraint = new BoundaryConstraint(
      this, "Intensity", std::numeric_limits<double>::epsilon(), true);
86
87
  addConstraint(IntensityConstraint);

88
89
90
  BoundaryConstraint *RadiusConstraint = new BoundaryConstraint(
      this, "Radius", std::numeric_limits<double>::epsilon(), true);
  addConstraint(RadiusConstraint);
91

92
93
94
  BoundaryConstraint *DecayConstraint = new BoundaryConstraint(
      this, "Decay", std::numeric_limits<double>::epsilon(), true);
  addConstraint(DecayConstraint);
95
96
}

97
98
99
100
101
102
103
void InelasticDiffRotDiscreteCircle::function1D(double *out,
                                                const double *xValues,
                                                const size_t nData) const {
  const double I = getParameter("Intensity");
  const double R = getParameter("Radius");
  const double rate = m_hbar / getParameter("Decay"); // micro-eV or mili-eV
  const int N = getAttribute("N").asInt();
104
  const double S = getParameter("Shift");
105

106
  double Q;
107
  if (getAttribute("Q").asDouble() == EMPTY_DBL()) {
108
    if (m_qValueCache.empty()) {
109
110
111
      throw std::runtime_error(
          "No Q attribute provided and cannot retrieve from worksapce.");
    }
112
113
114
115
116
117
    const int specIdx = getAttribute("WorkspaceIndex").asInt();
    Q = m_qValueCache[specIdx];

    g_log.debug() << "Get Q value for workspace index " << specIdx << ": " << Q
                  << std::endl;
  } else {
Dan Nixon's avatar
Dan Nixon committed
118
    Q = getAttribute("Q").asDouble();
119
120

    g_log.debug() << "Using Q attribute: " << Q << std::endl;
121
122
  }

123
124
125
126
127
  std::vector<double> sph(N);
  for (int k = 1; k < N; k++) {
    double x = 2 * Q * R * sin(M_PI * k / N);
    sph[k] =
        sin(x) / x; // spherical Besell function of order zero 'j0' is sin(x)/x
128
129
  }

130
131
  std::vector<double> ratel(N);
  for (int l = 1; l < N; l++) // l goes up to N-1
132
  {
133
    ratel[l] = rate * 4 * pow(sin(M_PI * l / N), 2); // notice that 0 < l/N < 1
134
135
  }

136
  for (size_t i = 0; i < nData; i++) {
137
    double w = xValues[i] - S;
138
    double S = 0.0;
139
    for (int l = 1; l < N; l++) // l goes up to N-1
140
    {
141
      double lorentzian = ratel[l] / (ratel[l] * ratel[l] + w * w);
142
      double al = 0.0;
143
      for (int k = 1; k < N; k++) // case k==N after the loop
144
145
      {
        double y = 2 * M_PI * l * k / N;
146
        al += cos(y) * sph[k];
147
148
149
150
151
      }
      al += 1; // limit for j0 when k==N, or x==0
      al /= N;
      S += al * lorentzian;
    }
152
    out[i] = I * S / M_PI;
153
154
155
  }
}

156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
/**
 * Handle seting fit workspace.
 *
 * Creates a list of Q values from each spectrum to be used with WorkspaceIndex
 * attribute.
 *
 * @param ws Pointer to workspace
 */
void InelasticDiffRotDiscreteCircle::setWorkspace(
    boost::shared_ptr<const API::Workspace> ws) {
  m_qValueCache.clear();

  auto workspace = boost::dynamic_pointer_cast<const MatrixWorkspace>(ws);
  if (!workspace)
    return;

  size_t numHist = workspace->getNumberHistograms();
  for (size_t idx = 0; idx < numHist; idx++) {
174
175
176
    IDetector_const_sptr det;
    try {
      det = workspace->getDetector(idx);
177
    } catch (Kernel::Exception::NotFoundError &) {
178
179
180
181
      m_qValueCache.clear();
      g_log.information("Cannot populate Q values from workspace");
      break;
    }
182

183
184
185
    try {
      double efixed = workspace->getEFixed(det);
      double usignTheta = workspace->detectorTwoTheta(det) / 2.0;
186

187
      double q = Mantid::Kernel::UnitConversion::run(usignTheta, efixed);
188

189
      m_qValueCache.push_back(q);
190
    } catch (std::runtime_error &) {
191
192
193
194
      m_qValueCache.clear();
      g_log.information("Cannot populate Q values from workspace");
      return;
    }
195
196
197
  }
}

198
199
200
201
/* Propagate the attribute to its member functions.
 * NOTE: we pass this->getAttribute(name) by reference, thus the same
 * object is shared by the composite function and its members.
 */
202
203
204
205
206
void DiffRotDiscreteCircle::trickleDownAttribute(const std::string &name) {
  for (size_t iFun = 0; iFun < nFunctions(); iFun++) {
    API::IFunction_sptr fun = getFunction(iFun);
    if (fun->hasAttribute(name))
      fun->setAttribute(name, this->getAttribute(name));
207
208
209
210
211
212
  }
}

/* Same as parent except we overwrite attributes of member functions
 * having the same name
 */
213
214
215
216
void DiffRotDiscreteCircle::declareAttribute(
    const std::string &name, const API::IFunction::Attribute &defaultValue) {
  API::ImmutableCompositeFunction::declareAttribute(name, defaultValue);
  trickleDownAttribute(name);
217
218
219
220
221
}

/* Same as parent except we overwrite attributes of member functions
 * having the same name
 */
222
223
224
225
void DiffRotDiscreteCircle::setAttribute(const std::string &name,
                                         const Attribute &att) {
  API::ImmutableCompositeFunction::setAttribute(name, att);
  trickleDownAttribute(name);
226
227
}

228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
// DiffRotDiscreteCircle::DiffRotDiscreteCircle()
void DiffRotDiscreteCircle::init() {
  m_elastic = boost::dynamic_pointer_cast<ElasticDiffRotDiscreteCircle>(
      API::FunctionFactory::Instance().createFunction(
          "ElasticDiffRotDiscreteCircle"));
  addFunction(m_elastic);
  m_inelastic = boost::dynamic_pointer_cast<InelasticDiffRotDiscreteCircle>(
      API::FunctionFactory::Instance().createFunction(
          "InelasticDiffRotDiscreteCircle"));
  addFunction(m_inelastic);

  setAttributeValue("NumDeriv", true);

  declareAttribute("Q", API::IFunction::Attribute(0.5));
  declareAttribute("N", API::IFunction::Attribute(3));

  // Set the aliases
  setAlias("f1.Intensity", "Intensity");
  setAlias("f1.Radius", "Radius");
  setAlias("f1.Decay", "Decay");
Dan Nixon's avatar
Dan Nixon committed
248
  setAlias("f1.Shift", "Shift");
249
250
251

  // Set the ties between Elastic and Inelastic parameters
  addDefaultTies("f0.Height=f1.Intensity,f0.Radius=f1.Radius");
252
  applyTies();
253
254
}

255
} // namespace Functions
256
257
} // namespace CurveFitting
} // namespace Mantid