GetDetectorOffsets.cpp 6.21 KB
Newer Older
1
#include "MantidAlgorithms/GetDetectorOffsets.h"
2
#include "MantidAPI/FileProperty.h"
3
#include "MantidAPI/FunctionFactory.h"
4
5
6
#include "MantidAPI/SpectraDetectorMap.h"
#include "MantidAPI/WorkspaceValidators.h"
#include "MantidDataObjects/OffsetsWorkspace.h"
7
#include <boost/math/special_functions/fpclassify.hpp>
8
9
#include <fstream>
#include <iomanip>
10
#include <ostream>
11
#include <sstream>
12

13

14
15
namespace Mantid
{
16
17
  namespace Algorithms
  {
18

19
20
    // Register the class into the algorithm factory
    DECLARE_ALGORITHM(GetDetectorOffsets)
21
22
23
24
    
    /// Sets documentation strings for this algorithm
    void GetDetectorOffsets::initDocs()
    {
25
26
      this->setWikiSummary("Creates an [[OffsetsWorkspace]] containing offsets for each detector. You can then save these to a .cal file using SaveCalFile.");
      this->setOptionalMessage("Creates an OffsetsWorkspace containing offsets for each detector. You can then save these to a .cal file using SaveCalFile.");
27
28
    }
    
29
30
    using namespace Kernel;
    using namespace API;
Peterson, Peter's avatar
Peterson, Peter committed
31
    using std::size_t;
32
    using namespace DataObjects;
33
34
35
36
37
38
39
40
41
42

    /// Constructor
    GetDetectorOffsets::GetDetectorOffsets() :
      API::Algorithm()
    {}

    /// Destructor
    GetDetectorOffsets::~GetDetectorOffsets()
    {}

43
44

    //-----------------------------------------------------------------------------------------
45
46
47
48
    /** Initialisation method. Declares properties to be used in algorithm.
     */
    void GetDetectorOffsets::init()
    {
49
50
51

      declareProperty(new WorkspaceProperty<>("InputWorkspace","",Direction::Input,
          new WorkspaceUnitValidator<>("dSpacing")),"A 2D workspace with X values of d-spacing");
52

53
54
55
56
57
58
59
60
61
      BoundedValidator<double> *mustBePositive = new BoundedValidator<double>();
      mustBePositive->setLower(0);

      declareProperty("Step",0.001, mustBePositive,
        "Step size used to bin d-spacing data");
      declareProperty("DReference",2.0, mustBePositive->clone(),
         "Center of reference peak in d-space");
      declareProperty("XMin",0.0, "Minimum of CrossCorrelation data to search for peak, usually negative");
      declareProperty("XMax",0.0, "Maximum of CrossCorrelation data to search for peak, usually positive");
62

63
64
      declareProperty(new FileProperty("GroupingFileName","", FileProperty::OptionalSave, ".cal"),
          "Optional: The name of the output CalFile to save the generated OffsetsWorkspace." );
65
66
67

      declareProperty(new WorkspaceProperty<OffsetsWorkspace>("OutputWorkspace","",Direction::Output),
          "An output workspace containing the offsets.");
68
69
    }

70
    //-----------------------------------------------------------------------------------------
71
72
73
74
75
76
77
78
79
    /** Executes the algorithm
     *
     *  @throw Exception::FileError If the grouping file cannot be opened or read successfully
     */
    void GetDetectorOffsets::exec()
    {
      inputW=getProperty("InputWorkspace");
      Xmin=getProperty("XMin");
      Xmax=getProperty("XMax");
80
81
      if (Xmin>=Xmax)
        throw std::runtime_error("Must specify Xmin<Xmax");
82
83
84
      dreference=getProperty("DReference");
      step=getProperty("Step");
      nspec=inputW->getNumberHistograms();
85
86
87
88
89
      // Create the output OffsetsWorkspace
      OffsetsWorkspace_sptr outputW(new OffsetsWorkspace(inputW->getInstrument()));

      // Fit all the spectra with a gaussian
      Progress prog(this, 0, 1.0, nspec);
90
      PARALLEL_FOR1(inputW)
91
92
      for (int i=0;i<nspec;++i)
      {
93
        PARALLEL_START_INTERUPT_REGION
94
95
        double offset=fitSpectra(i);
        int detID = inputW->getDetector(i)->getID();
96
97
98
99
100

        PARALLEL_CRITICAL(GetDetectorOffsets_setValue)
        { // Most of the exec time is in FitSpectra, so this critical block should not be a problem.
          outputW->setValue(detID, offset);
        }
101
        prog.report();
102
        PARALLEL_END_INTERUPT_REGION
103
      }
104
      PARALLEL_CHECK_INTERUPT_REGION
105
106
107

      // Return the output
      setProperty("OutputWorkspace",outputW);
108
109
110
111
112
113
114
115
116
117
118
119

      // Also save to .cal file, if requested
      std::string filename=getProperty("GroupingFileName");
      if (!filename.empty())
      {
        progress(0.9, "Saving .cal file");
        IAlgorithm_sptr childAlg = createSubAlgorithm("SaveCalFile");
        childAlg->setProperty("OffsetsWorkspace", outputW);
        childAlg->setPropertyValue("Filename", filename);
        childAlg->executeAsSubAlg();
      }

120
121
    }

122
123

    //-----------------------------------------------------------------------------------------
124
   /** Calls Gaussian1D as a child algorithm to fit the offset peak in a spectrum
125
126
    *
    *  @param s :: The Workspace Index to fit
127
128
    *  @return The calculated offset value
    */
Peterson, Peter's avatar
Peterson, Peter committed
129
    double GetDetectorOffsets::fitSpectra(const int64_t s)
130
131
132
133
134
135
    {
      // Find point of peak centre
      const MantidVec & yValues = inputW->readY(s);
      MantidVec::const_iterator it = std::max_element(yValues.begin(), yValues.end());
      const double peakHeight = *it; 
      const double peakLoc = inputW->readX(s)[it - yValues.begin()];
136
      // Return offset of 0 if peak of Cross Correlation is nan (Happens when spectra is zero)
137
      if ( peakHeight < 0.75 || boost::math::isnan(peakHeight) ) return (0.);
138
139
140
141

      IAlgorithm_sptr fit_alg;
      try
      {
142
        //set the subalgorithm no to log as this will be run once per spectra
143
        fit_alg = createSubAlgorithm("Fit",-1,-1,false);
144
145
      } catch (Exception::NotFoundError&)
      {
146
        g_log.error("Can't locate Fit algorithm");
147
        throw ;
148
149
      }
      fit_alg->setProperty("InputWorkspace",inputW);
150
      fit_alg->setProperty<int>("WorkspaceIndex",s);
151
152
153
154
      fit_alg->setProperty("StartX",Xmin);
      fit_alg->setProperty("EndX",Xmax);
      fit_alg->setProperty("MaxIterations",100);

155
156
157
158
159
      std::ostringstream fun_str;
      fun_str << "name=LinearBackground;name=Gaussian,Height="<<peakHeight<<",";
      fun_str << "PeakCentre="<<peakLoc<<",Sigma=10.0";

      fit_alg->setProperty("Function",fun_str.str());
160
      fit_alg->executeAsSubAlg();
161
      std::string fitStatus = fit_alg->getProperty("OutputStatus");
162
      if ( fitStatus.compare("success") ) return (0.);
163

Roman Tolchenov's avatar
Roman Tolchenov committed
164
165
      std::vector<double> params = fit_alg->getProperty("Parameters");
      const double offset = params[3]; // f1.PeakCentre
Peterson, Peter's avatar
Peterson, Peter committed
166
      return (-1.*offset*step/(dreference+offset*step));
167
168
169
170
171
    }



  } // namespace Algorithm
172
} // namespace Mantid