CrystalFieldSpectrum.cpp 5.62 KB
Newer Older
1
#include "MantidCurveFitting/Functions/CrystalFieldSpectrum.h"
2
#include "MantidCurveFitting/Functions/CrystalFieldPeaks.h"
3
#include "MantidCurveFitting/Functions/CrystalFieldPeakUtils.h"
4
5

#include "MantidAPI/CompositeFunction.h"
6
#include "MantidAPI/IConstraint.h"
7
#include "MantidAPI/IPeakFunction.h"
8
#include "MantidAPI/FunctionFactory.h"
9
#include "MantidAPI/ParameterTie.h"
10
#include "MantidCurveFitting/Constraints/BoundaryConstraint.h"
11
#include "MantidKernel/Strings.h"
12
13
14

#include <algorithm>
#include <iostream>
15
16
17
18
19
20
21
22
23
24
25
26

namespace Mantid {
namespace CurveFitting {
namespace Functions {

using namespace CurveFitting;
using namespace Kernel;
using namespace API;

DECLARE_FUNCTION(CrystalFieldSpectrum)

/// Constructor
27
CrystalFieldSpectrum::CrystalFieldSpectrum()
28
    : FunctionGenerator(IFunction_sptr(new CrystalFieldPeaks)), m_nPeaks(0) {
29
  declareAttribute("PeakShape", Attribute("Lorentzian"));
30
  declareAttribute("FWHM", Attribute(0.0));
31
  std::vector<double> vec;
32
33
34
  declareAttribute("FWHMX", Attribute(vec));
  declareAttribute("FWHMY", Attribute(vec));
  declareAttribute("FWHMVariation", Attribute(0.1));
35
  declareAttribute("NPeaks", Attribute(0));
36
  declareAttribute("FixAllPeaks", Attribute(false));
37
38
}

39
40
/// Uses m_crystalField to calculate peak centres and intensities
/// then populates m_spectrum with peaks of type given in PeakShape attribute.
41
void CrystalFieldSpectrum::buildTargetFunction() const {
42
  m_dirty = false;
43
44
  auto spectrum = new CompositeFunction;
  m_target.reset(spectrum);
45
  m_target->setAttribute("NumDeriv", this->getAttribute("NumDeriv"));
46

47
48
  FunctionDomainGeneral domain;
  FunctionValues values;
49
  m_source->function(domain, values);
50
51
52
53
54
55

  if (values.size() == 0) {
    return;
  }

  if (values.size() % 2 != 0) {
56
    throw std::runtime_error(
57
        "CrystalFieldPeaks returned odd number of values.");
58
59
  }

60
61
  auto xVec = getAttribute("FWHMX").asVector();
  auto yVec = getAttribute("FWHMY").asVector();
62
  auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
63

64
65
66
67
  auto peakShape = getAttribute("PeakShape").asString();
  auto defaultFWHM = getAttribute("FWHM").asDouble();
  size_t nRequiredPeaks = getAttribute("NPeaks").asInt();
  bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();
68
  m_nPeaks = CrystalFieldUtils::buildSpectrumFunction(
69
      *spectrum, peakShape, values, xVec, yVec, fwhmVariation, defaultFWHM,
70
      nRequiredPeaks, fixAllPeaks);
71
  storeReadOnlyAttribute("NPeaks", Attribute(static_cast<int>(m_nPeaks)));
72
}
73

74
/// Update m_spectrum function.
75
76
77
void CrystalFieldSpectrum::updateTargetFunction() const {
  if (!m_target) {
    buildTargetFunction();
78
79
80
    return;
  }
  m_dirty = false;
81
  auto peakShape = getAttribute("PeakShape").asString();
82
83
  auto xVec = getAttribute("FWHMX").asVector();
  auto yVec = getAttribute("FWHMY").asVector();
84
  auto fwhmVariation = getAttribute("FWHMVariation").asDouble();
85
86
  auto defaultFWHM = getAttribute("FWHM").asDouble();
  bool fixAllPeaks = getAttribute("FixAllPeaks").asBool();
87
88
  FunctionDomainGeneral domain;
  FunctionValues values;
89
  m_source->function(domain, values);
90
  m_target->setAttribute("NumDeriv", this->getAttribute("NumDeriv"));
91
  auto &spectrum = dynamic_cast<CompositeFunction &>(*m_target);
92
  m_nPeaks = CrystalFieldUtils::calculateNPeaks(values);
93
94
95
  CrystalFieldUtils::updateSpectrumFunction(
      spectrum, peakShape, values, m_nPeaks, 0, xVec, yVec, fwhmVariation,
      defaultFWHM, fixAllPeaks);
96
  storeReadOnlyAttribute("NPeaks", Attribute(static_cast<int>(m_nPeaks)));
97
}
98

99
100
101
102
103
104
105
106
/// Custom string conversion method
std::string CrystalFieldSpectrum::asString() const {
  std::ostringstream ostr;
  ostr << "name=" << this->name();
  // Print the attributes
  std::vector<std::string> attr = this->getAttributeNames();
  for (const auto &attName : attr) {
    std::string attValue = this->getAttribute(attName).value();
107
    if (!attValue.empty() && attValue != "\"\"" && attValue != "()") {
108
109
110
      ostr << ',' << attName << '=' << attValue;
    }
  }
111
  std::vector<std::string> ties;
112
  // Print own parameters
113
  for (size_t i = 0; i < m_nOwnParams; i++) {
114
115
116
117
118
119
    std::ostringstream paramOut;
    paramOut << parameterName(i) << '=' << getParameter(i);
    if (isActive(i)) {
      ostr << ',' << paramOut.str();
    } else if (isFixed(i)) {
      ties.push_back(paramOut.str());
120
121
122
    }
  }

123
124
125
126
  // Print parameters of the important peaks only
  const CompositeFunction &spectrum =
      dynamic_cast<const CompositeFunction &>(*m_target);
  for (size_t ip = 0; ip < m_nPeaks; ++ip) {
127
    const auto &peak = dynamic_cast<IPeakFunction &>(*spectrum.getFunction(ip));
128
129
130
131
132
133
134
135
136
137
138
139
    // Print peak's atributes
    auto attr = peak.getAttributeNames();
    for (const auto &attName : attr) {
      std::string attValue = peak.getAttribute(attName).value();
      if (!attValue.empty() && attValue != "\"\"") {
        ostr << ",f" << ip << "." << attName << '=' << attValue;
      }
    }
    // Print peak's parameters
    for (size_t i = 0; i < peak.nParams(); i++) {
      const ParameterTie *tie = peak.getTie(i);
      if (!tie || !tie->isDefault()) {
140
141
        ostr << ",f" << ip << "." << peak.parameterName(i) << '='
             << peak.getParameter(i);
142
      }
143
144
145
    }
  } // for peaks

146
147
  // collect non-default constraints
  std::string constraints = writeConstraints();
148
149
  // print constraints
  if (!constraints.empty()) {
Roman Tolchenov's avatar
Roman Tolchenov committed
150
    ostr << ",constraints=(" << constraints << ")";
151
152
  }

153
154
155
156
157
  // collect the non-default ties
  auto tiesString = writeTies();
  if (!tiesString.empty()) {
    ties.push_back(tiesString);
  }
158
159
160
161
  // print the ties
  if (!ties.empty()) {
    ostr << ",ties=(" << Kernel::Strings::join(ties.begin(), ties.end(), ",")
         << ")";
162
163
164
165
166
  }

  return ostr.str();
}

167
168
169
} // namespace Functions
} // namespace CurveFitting
} // namespace Mantid