NormaliseVanadium.cpp 4.94 KB
Newer Older
1
2
3
4
//----------------------------------------------------------------------
// Includes
//----------------------------------------------------------------------
#include "MantidCrystal/NormaliseVanadium.h"
5
#include "MantidAPI/Axis.h"
6
#include "MantidAPI/InstrumentValidator.h"
7
#include "MantidGeometry/Instrument.h"
8
#include "MantidKernel/BoundedValidator.h"
9
10
11
12
13
14
15
#include "MantidKernel/UnitFactory.h"
#include "MantidKernel/Fast_Exponential.h"
#include "MantidKernel/VectorHelper.h"

/*  Following A.J.Schultz's anvred, scaling the vanadium spectra:
 */

16
17
namespace Mantid {
namespace Crystal {
18
19
20
21
22
23
24
25
26

// Register the class into the algorithm factory
DECLARE_ALGORITHM(NormaliseVanadium)

using namespace Kernel;
using namespace Geometry;
using namespace API;
using namespace DataObjects;

27
NormaliseVanadium::NormaliseVanadium() : API::Algorithm() {}
28

29
void NormaliseVanadium::init() {
30
  // The input workspace must have an instrument and units of wavelength
31
  auto wsValidator = boost::make_shared<InstrumentValidator>();
32

33
34
35
36
37
38
39
  declareProperty(new WorkspaceProperty<>("InputWorkspace", "",
                                          Direction::Input, wsValidator),
                  "The X values for the input workspace must be in units of "
                  "wavelength or TOF");
  declareProperty(
      new WorkspaceProperty<>("OutputWorkspace", "", Direction::Output),
      "Output workspace name");
40

41
  auto mustBePositive = boost::make_shared<BoundedValidator<double>>();
42
43
  mustBePositive->setLower(0.0);
  declareProperty("Wavelength", 1.0, mustBePositive,
44
                  "Normalizes spectra to this wavelength");
45
46
}

47
void NormaliseVanadium::exec() {
48
49
50
51
52
53
54
  // 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

55
56
  MatrixWorkspace_sptr correctionFactors =
      WorkspaceFactory::Instance().create(m_inputWS);
57

58
59
  const int64_t numHists =
      static_cast<int64_t>(m_inputWS->getNumberHistograms());
60
61
62
63
64
65
  const int64_t specSize = static_cast<int64_t>(m_inputWS->blocksize());

  const bool isHist = m_inputWS->isHistogramData();

  // If sample not at origin, shift cached positions.
  const V3D samplePos = m_inputWS->getInstrument()->getSample()->getPos();
66
  const V3D pos = m_inputWS->getInstrument()->getSource()->getPos() - samplePos;
67
68
  double L1 = pos.norm();

69
  Progress prog(this, 0.0, 1.0, numHists);
70
  // Loop over the spectra
71
72
73
  PARALLEL_FOR2(m_inputWS, correctionFactors)
  for (int64_t i = 0; i < int64_t(numHists); ++i) {
    //    PARALLEL_START_INTERUPT_REGION //FIXME: Restore
74
75

    // Get a reference to the Y's in the output WS for storing the factors
76
77
    MantidVec &Y = correctionFactors->dataY(i);
    MantidVec &E = correctionFactors->dataE(i);
78
79

    // Copy over bin boundaries
80
    const ISpectrum *inSpec = m_inputWS->getSpectrum(i);
81
    inSpec->lockData(); // for MRU-related thread safety
82

83
    const MantidVec &Xin = inSpec->readX();
84
    correctionFactors->dataX(i) = Xin;
85
86
    const MantidVec &Yin = inSpec->readY();
    const MantidVec &Ein = inSpec->readE();
87
88
89

    // Get detector position
    IDetector_const_sptr det;
90
    try {
91
      det = m_inputWS->getDetector(i);
92
93
94
95
96
    } catch (Exception::NotFoundError &) {
      // Catch if no detector. Next line tests whether this happened - test
      // placed
      // outside here because Mac Intel compiler doesn't like 'continue' in a
      // catch
97
98
99
      // in an openmp block.
    }
    // If no detector found, skip onto the next spectrum
100
101
    if (!det)
      continue;
102
103
104
105
106

    // This is the scattered beam direction
    Instrument_const_sptr inst = m_inputWS->getInstrument();
    V3D dir = det->getPos() - samplePos;
    double L2 = dir.norm();
107
108
109
    // Two-theta = polar angle = scattering angle = between +Z vector and the
    // scattered beam
    double scattering = dir.angle(V3D(0.0, 0.0, 1.0));
110
111
112
113
114
115
116
117
118

    Mantid::Kernel::Units::Wavelength wl;
    std::vector<double> timeflight;

    // Loop through the bins in the current spectrum
    double lambp = 0;
    double lambm = 0;
    double normp = 0;
    double normm = 0;
119
    for (int64_t j = 0; j < specSize; j++) {
120
121
122
123
      timeflight.push_back((isHist ? (0.5 * (Xin[j] + Xin[j + 1])) : Xin[j]));
      if (unitStr.compare("TOF") == 0)
        wl.fromTOF(timeflight, timeflight, L1, L2, scattering, 0, 0, 0);
      const double lambda = timeflight[0];
124
      if (lambda > lambdanorm) {
125
126
127
128
129
130
131
132
        lambp = lambda;
        normp = Yin[j];
        break;
      }
      lambm = lambda;
      normm = Yin[j];
      timeflight.clear();
    }
133
134
135
136
137
    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;
138
139
    }

140
    inSpec->unlockData();
141
142
143

    prog.report();

144
    //    PARALLEL_END_INTERUPT_REGION
145
  }
146
  //  PARALLEL_CHECK_INTERUPT_REGION
147
148
149
150
151
152

  setProperty("OutputWorkspace", correctionFactors);
}

} // namespace Crystal
} // namespace Mantid