LoadPreNexusMonitors.cpp 9.28 KB
Newer Older
1
#include "MantidDataHandling/LoadPreNexusMonitors.h"
2
3

#include "MantidAPI/Axis.h"
4
#include "MantidAPI/FileProperty.h"
5
#include "MantidAPI/MatrixWorkspace.h"
6
#include "MantidAPI/WorkspaceFactory.h"
7
#include "MantidKernel/BinaryFile.h"
8
9
10
#include "MantidKernel/ConfigService.h"
#include "MantidKernel/UnitFactory.h"

11
12
13
14
#include <boost/lexical_cast.hpp>
#include <boost/shared_array.hpp>

#include <Poco/DOM/AutoPtr.h>
Campbell, Stuart's avatar
Campbell, Stuart committed
15
16
17
18
19
20
#include <Poco/DOM/DOMParser.h>
#include <Poco/DOM/Document.h>
#include <Poco/DOM/Element.h>
#include <Poco/DOM/NodeIterator.h>
#include <Poco/DOM/NodeFilter.h>
#include <Poco/DOM/NodeList.h>
21
#include <Poco/Path.h>
Campbell, Stuart's avatar
Campbell, Stuart committed
22
23
#include <Poco/SAX/InputSource.h>

24
25
26
27
28
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iterator>

29
30
namespace Mantid {
namespace DataHandling {
31
32

// Register the algorithm into the AlgorithmFactory
33
DECLARE_ALGORITHM(LoadPreNexusMonitors)
34
35
36
37
38
39
40
41
42
43
44

using namespace Mantid::Kernel;
using namespace Mantid::API;

// Some constants for property names
static const std::string RUNINFO_FILENAME("RunInfoFilename");
static const std::string WORKSPACE_OUT("OutputWorkspace");

// A reference to the logger is provided by the base class, it is called g_log.
// It is used to print out information, warning and error messages

45
46
47
48
LoadPreNexusMonitors::LoadPreNexusMonitors()
    : Mantid::API::Algorithm(), nMonitors(0),
      instrument_loaded_correctly(false) {}

49
void LoadPreNexusMonitors::init() {
50
  // Filename for the runinfo file.
51
52
53
54
  declareProperty(new FileProperty(RUNINFO_FILENAME, "", FileProperty::Load,
                                   "_runinfo.xml"),
                  "The filename of the runinfo file for a particular run. "
                  "Allowed Values are: _runinfo.xml");
55
56

  // The output workspace
57
58
  declareProperty(new WorkspaceProperty<MatrixWorkspace>(WORKSPACE_OUT, "",
                                                         Direction::Output),
59
                  "The workspace to load the monitors into.");
60
61
62
63
64

  // Make sure things are initialised.
  nMonitors = 0;
}

65
void LoadPreNexusMonitors::exec() {
66
  // time of flight channel parameters
67
68
69
  double tmin = 0.0;
  double tstep = 0.0;
  int tchannels = 0;
70
71
72
  std::string instrumentName;

  // Vectors to store monitor parameters
73
74
  std::vector<std::string> monitorFilenames;
  std::vector<int> monitorIDs;
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

  // Get the Runinfo filename from the property
  std::string runinfo_filename = this->getPropertyValue(RUNINFO_FILENAME);

  // TODO: Extract the directory that the runinfo file is in.
  // Create a Poco Path object for runinfo filename
  Poco::Path runinfoPath(runinfo_filename, Poco::Path::PATH_GUESS);
  // Now lets get the directory
  Poco::Path dirPath(runinfoPath.parent());

  this->g_log.information("Monitor File Dir: " + dirPath.toString());

  // Some XML parsing magic...
  std::ifstream in(runinfo_filename.c_str());
  Poco::XML::InputSource src(in);

  Poco::XML::DOMParser parser;
  Poco::AutoPtr<Poco::XML::Document> pDoc = parser.parse(&src);

  Poco::XML::NodeIterator it(pDoc, Poco::XML::NodeFilter::SHOW_ELEMENT);
95
96
  Poco::XML::Node *pNode = it.nextNode();
  while (pNode) {
97
98

    // Get the beamline name.
99
100
101
    if (pNode->nodeName() == "RunInfo") {
      Poco::XML::Element *pRunInfoElement =
          static_cast<Poco::XML::Element *>(pNode);
102
103
104
105
      instrumentName = pRunInfoElement->getAttribute("instrument");
    }

    // Look for the section that contains the binning parameters.
106
107
    // This will end up with using the binning for the last monitor that is
    // found
108
    // TODO: Modify to store (and use) the TCBs for each monitor separately.
109
    if (pNode->nodeName() == "BeamMonitorInfo") {
110
111
112
      // Increment the number of monitors we've found
      ++nMonitors;

113
      Poco::XML::Element *pE = static_cast<Poco::XML::Element *>(pNode);
114
115
      g_log.debug() << "Beam Monitor " << pE->getAttribute("id") << std::endl;
      g_log.debug() << "\tname: " << pE->getAttribute("name") << std::endl;
116
117
      g_log.debug() << "\tdescription: " << pE->getAttribute("description")
                    << std::endl;
118

119
      // Now lets get the tof binning settings
120
121
122
123
      Poco::XML::Element *pTimeChannels =
          pE->getChildElement("NumTimeChannels");
      tmin =
          boost::lexical_cast<double>(pTimeChannels->getAttribute("startbin"));
124
      tstep = boost::lexical_cast<double>(pTimeChannels->getAttribute("width"));
125
126
127
128
    }

    // Look for the 'DataList' node to get the monitor dims.
    // TODO: Again we will only use the mast monitor value.
129
    if (pNode->nodeName() == "DataList") {
130
      // Get a list of the child elements
131
132
133
      Poco::AutoPtr<Poco::XML::NodeList> pDataListChildren =
          pNode->childNodes();
      for (unsigned long i = 0; i < pDataListChildren->length(); ++i) {
134
        // We only care about monitors
135
136
137
138
139
140
        if (pDataListChildren->item(i)->nodeName() == "monitor") {
          Poco::XML::Element *element =
              static_cast<Poco::XML::Element *>(pDataListChildren->item(i));
          monitorIDs.push_back(
              boost::lexical_cast<int>(element->getAttribute("id")));
          monitorFilenames.push_back(element->getAttribute("name"));
141
142
143
144
145
        }
      }
    }

    // Get the size of the files
146
    if (pNode->nodeName() == "FileFormats") {
147
      // Get a list of the child elements
148
149
150
      Poco::AutoPtr<Poco::XML::NodeList> pDataListChildren =
          pNode->childNodes();
      for (unsigned long i = 0; i < pDataListChildren->length(); ++i) {
151
        // We only care about monitors
152
        if (pDataListChildren->item(i)->nodeName() == "monitor") {
153
          std::string dims =
154
155
              static_cast<Poco::XML::Element *>(pDataListChildren->item(i))
                  ->getAttribute("dims");
156
          tchannels = boost::lexical_cast<int>(dims);
157
158
159
160
161
162
163
        }
      }
    }

    pNode = it.nextNode();
  }

164
165
  g_log.information() << "Found " << nMonitors << " beam monitors."
                      << std::endl;
166

167
168
  g_log.information() << "Number of Time Channels = " << tchannels << std::endl;

169
  // Now lets create the time of flight array.
170
171
  const int numberTimeBins = tchannels + 1;
  MantidVec time_bins(numberTimeBins);
172
173
  for (int i = 0; i < numberTimeBins; ++i) {
    time_bins[i] = tmin + (i)*tstep;
174
175
176
  }

  // Create the new workspace
177
178
  MatrixWorkspace_sptr localWorkspace = WorkspaceFactory::Instance().create(
      "Workspace2D", nMonitors, numberTimeBins, tchannels);
179

180
  for (int i = 0; i < nMonitors; i++) {
181
182
183
    // Now lets actually read the monitor files..
    Poco::Path pMonitorFilename(dirPath, monitorFilenames[i]);

184
185
    g_log.debug() << "Loading monitor file :" << pMonitorFilename.toString()
                  << std::endl;
186

187
    Kernel::BinaryFile<uint32_t> monitorFile(pMonitorFilename.toString());
188
189
    // temp buffer for file reading
    std::vector<uint32_t> buffer = monitorFile.loadAllIntoVector();
190
191

    MantidVec intensity(buffer.begin(), buffer.end());
192
    // Copy the same data into the error array
193
    MantidVec error(buffer.begin(), buffer.end());
194
    // Now take the sqrt()
195
196
    std::transform(error.begin(), error.end(), error.begin(),
                   (double (*)(double))sqrt);
197
198
199
200
201

    localWorkspace->dataX(i) = time_bins;
    localWorkspace->dataY(i) = intensity;
    localWorkspace->dataE(i) = error;
    // Just have spectrum number be the same as the monitor number but -ve.
202
    ISpectrum *spectrum = localWorkspace->getSpectrum(i);
203
204
    spectrum->setSpectrumNo(monitorIDs[i]);
    spectrum->setDetectorID(-monitorIDs[i]);
205
206
  }

207
  g_log.debug() << "Setting axis zero to TOF" << std::endl;
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222

  // Set the units
  localWorkspace->getAxis(0)->unit() = UnitFactory::Instance().create("TOF");
  localWorkspace->setYUnit("Counts");
  // TODO localWorkspace->setTitle(title);

  // Actually load the instrument
  this->runLoadInstrument(instrumentName, localWorkspace);

  // Set the property
  setProperty("OutputWorkspace", localWorkspace);
}

//-----------------------------------------------------------------------------
/** Load the instrument geometry File
223
 *  @param instrument :: instrument name.
224
225
 *  @param localWorkspace :: MatrixWorkspace in which to put the instrument
 * geometry
226
 */
227
228
void LoadPreNexusMonitors::runLoadInstrument(
    const std::string &instrument, MatrixWorkspace_sptr localWorkspace) {
229

230
  IAlgorithm_sptr loadInst = createChildAlgorithm("LoadInstrument");
231

232
  // Now execute the Child Algorithm. Catch and log any error, but don't stop.
233
  bool executionSuccessful(true);
234
  try {
235
    loadInst->setPropertyValue("InstrumentName", instrument);
236
237
238
    loadInst->setProperty<MatrixWorkspace_sptr>("Workspace", localWorkspace);
    loadInst->setProperty("RewriteSpectraMap",
                          false); // We have a custom mapping
239
240
    loadInst->execute();

241
242
    // Populate the instrument parameters in this workspace - this works around
    // a bug
243
    localWorkspace->populateInstrumentParameters();
244
245
246
  } catch (std::invalid_argument &e) {
    g_log.information()
        << "Invalid argument to LoadInstrument Child Algorithm : " << e.what()
247
248
        << std::endl;
    executionSuccessful = false;
249
250
251
252
  } catch (std::runtime_error &e) {
    g_log.information()
        << "Unable to successfully run LoadInstrument Child Algorithm : "
        << e.what() << std::endl;
253
254
255
256
    executionSuccessful = false;
  }

  // If loading instrument definition file fails
257
  if (!executionSuccessful) {
258
    g_log.error() << "Error loading Instrument definition file\n";
259
  } else {
260
261
262
263
264
265
    this->instrument_loaded_correctly = true;
  }
}

} // namespace DataHandling
} // namespace Mantid