LoadIsawDetCal.cpp 16.7 KB
Newer Older
1
#include "MantidDataHandling/LoadIsawDetCal.h"
2

Owen Arnold's avatar
Owen Arnold committed
3
#include "MantidGeometry/Instrument/ComponentInfo.h"
4
#include "MantidAPI/FileProperty.h"
5
#include "MantidAPI/InstrumentValidator.h"
6
#include "MantidAPI/MultipleFileProperty.h"
7
#include "MantidAPI/Run.h"
8

9
#include "MantidDataObjects/EventList.h"
Peterson, Peter's avatar
Peterson, Peter committed
10
#include "MantidDataObjects/EventWorkspace.h"
11
#include "MantidDataObjects/PeaksWorkspace.h"
Peterson, Peter's avatar
Peterson, Peter committed
12
#include "MantidDataObjects/Workspace2D.h"
13

Peterson, Peter's avatar
Peterson, Peter committed
14
#include "MantidGeometry/Instrument.h"
Peterson, Peter's avatar
Peterson, Peter committed
15
16
#include "MantidGeometry/Instrument/ObjCompAssembly.h"
#include "MantidGeometry/Instrument/RectangularDetector.h"
17

18
#include "MantidKernel/Strings.h"
19
#include "MantidKernel/V3D.h"
20

21
#include <algorithm>
22
#include <boost/algorithm/string/trim.hpp>
23
#include <fstream>
24
#include <iostream>
25
#include <numeric>
26
#include <sstream>
27

28
29
namespace Mantid {
namespace DataHandling {
30

31
32
// Register the class into the algorithm factory
DECLARE_ALGORITHM(LoadIsawDetCal)
33

34
35
36
37
using namespace Kernel;
using namespace API;
using namespace Geometry;
using namespace DataObjects;
38

39
40
41
/** Initialisation method
*/
void LoadIsawDetCal::init() {
Peterson, Peter's avatar
Peterson, Peter committed
42
  declareProperty(Kernel::make_unique<WorkspaceProperty<Workspace>>(
43
44
45
                      "InputWorkspace", "", Direction::InOut,
                      boost::make_shared<InstrumentValidator>()),
                  "The workspace containing the geometry to be calibrated.");
46

47
  const auto exts = std::vector<std::string>({".DetCal"});
Peterson, Peter's avatar
Peterson, Peter committed
48
49
50
51
  declareProperty(
      Kernel::make_unique<API::MultipleFileProperty>("Filename", exts),
      "The input filename of the ISAW DetCal file (Two files "
      "allowed for SNAP) ");
52

53
54
55
56
57
58
  declareProperty(
      Kernel::make_unique<API::FileProperty>(
          "Filename2", "", API::FileProperty::OptionalLoad, ".DetCal"),
      "The input filename of the second ISAW DetCal file (West "
      "banks for SNAP) ");

59
60
61
  declareProperty("TimeOffset", 0.0, "Time Offset", Direction::Output);
}

62
63
namespace {
const constexpr double DegreesPerRadian = 180.0 / M_PI;
64

Peterson, Peter's avatar
Peterson, Peter committed
65
std::string getBankName(const std::string &bankPart, const int idnum) {
66
67
68
69
70
71
  if (bankPart == "WISHpanel" && idnum < 10) {
    return bankPart + "0" + std::to_string(idnum);
  } else {
    return bankPart + std::to_string(idnum);
  }
}
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

std::string getInstName(API::Workspace_const_sptr wksp) {
  MatrixWorkspace_const_sptr matrixWksp =
      boost::dynamic_pointer_cast<const MatrixWorkspace>(wksp);
  if (matrixWksp) {
    return matrixWksp->getInstrument()->getName();
  }

  PeaksWorkspace_const_sptr peaksWksp =
      boost::dynamic_pointer_cast<const PeaksWorkspace>(wksp);
  if (peaksWksp) {
    return peaksWksp->getInstrument()->getName();
  }

  throw std::runtime_error("Failed to determine instrument name");
}
}

std::map<std::string, std::string> LoadIsawDetCal::validateInputs() {
  std::map<std::string, std::string> result;

  // two detcal files is only valid for snap
  std::vector<std::string> filenames = getFilenames();
95
  if (filenames.empty()) {
96
97
98
99
100
    result["Filename"] = "Must supply .detcal file";
  } else if (filenames.size() == 2) {
    Workspace_const_sptr wksp = getProperty("InputWorkspace");
    const auto instname = getInstName(wksp);

101
    if (instname != "SNAP") {
102
103
104
      result["Filename"] = "Two files is only valid for SNAP";
    }
  } else if (filenames.size() > 2) {
105
    result["Filename"] = "Supply at most two .detcal files";
106
107
108
  }

  return result;
109
110
}

111
112
113
114
115
116
/** Executes the algorithm
*
*  @throw runtime_error Thrown if algorithm cannot execute
*/
void LoadIsawDetCal::exec() {
  // Get the input workspace
117
  Workspace_sptr ws = getProperty("InputWorkspace");
118
119
  MatrixWorkspace_sptr inputW =
      boost::dynamic_pointer_cast<MatrixWorkspace>(ws);
120
  PeaksWorkspace_sptr inputP = boost::dynamic_pointer_cast<PeaksWorkspace>(ws);
121

122
  Instrument_sptr inst = getCheckInst(ws);
123

124
125
  std::string instname = inst->getName();

126
  const auto filenames = getFilenames();
127
128

  // Output summary to log file
Hahn, Steven's avatar
Hahn, Steven committed
129
  int count, id, nrows, ncols;
130
131
  double width, height, depth, detd, x, y, z, base_x, base_y, base_z, up_x,
      up_y, up_z;
132
  std::ifstream input(filenames[0].c_str(), std::ios_base::in);
133
134
135
  std::string line;
  std::string detname;
  // Build a list of Rectangular Detectors
136
  std::vector<boost::shared_ptr<RectangularDetector>> detList;
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
  for (int i = 0; i < inst->nelements(); i++) {
    boost::shared_ptr<RectangularDetector> det;
    boost::shared_ptr<ICompAssembly> assem;
    boost::shared_ptr<ICompAssembly> assem2;

    det = boost::dynamic_pointer_cast<RectangularDetector>((*inst)[i]);
    if (det) {
      detList.push_back(det);
    } else {
      // Also, look in the first sub-level for RectangularDetectors (e.g. PG3).
      // We are not doing a full recursive search since that will be very long
      // for lots of pixels.
      assem = boost::dynamic_pointer_cast<ICompAssembly>((*inst)[i]);
      if (assem) {
        for (int j = 0; j < assem->nelements(); j++) {
          det = boost::dynamic_pointer_cast<RectangularDetector>((*assem)[j]);
          if (det) {
            detList.push_back(det);

          } else {
            // Also, look in the second sub-level for RectangularDetectors (e.g.
            // PG3).
            // We are not doing a full recursive search since that will be very
            // long for lots of pixels.
            assem2 = boost::dynamic_pointer_cast<ICompAssembly>((*assem)[j]);
            if (assem2) {
              for (int k = 0; k < assem2->nelements(); k++) {
                det = boost::dynamic_pointer_cast<RectangularDetector>(
                    (*assem2)[k]);
                if (det) {
                  detList.push_back(det);
168
169
170
                }
              }
            }
171
172
173
174
          }
        }
      }
    }
175
  }
176
  std::unordered_set<int> uniqueBanks; // for CORELLI and WISH
177
  std::string bankPart = "bank";
178
  if (instname == "WISH")
179
180
    bankPart = "WISHpanel";
  if (detList.empty()) {
181
182
183
184
    // Get all children
    std::vector<IComponent_const_sptr> comps;
    inst->getChildren(comps, true);

185
186
    for (auto &comp : comps) {
      std::string bankName = comp->getName();
187
      boost::trim(bankName);
188
      boost::erase_all(bankName, bankPart);
189
190
      int bank = 0;
      Strings::convert(bankName, bank);
191
192
      if (bank == 0)
        continue;
193
194
195
196
      // Track unique bank numbers
      uniqueBanks.insert(bank);
    }
  }
197

198
  auto expInfoWS = boost::dynamic_pointer_cast<ExperimentInfo>(ws);
199
  auto &componentInfo = expInfoWS->mutableComponentInfo();
200
  std::vector<ComponentScaling> rectangularDetectorScalings;
201

202
203
204
205
206
207
  while (std::getline(input, line)) {
    if (line[0] == '7') {
      double mL1, mT0;
      std::stringstream(line) >> count >> mL1 >> mT0;
      setProperty("TimeOffset", mT0);
      // Convert from cm to m
208
      if (instname == "WISH")
209
        center(0.0, 0.0, -mL1, "undulator", ws, componentInfo);
210
      else
211
        center(0.0, 0.0, -mL1, "moderator", ws, componentInfo);
212
      // mT0 and time of flight are both in microsec
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
      if (mT0 != 0.0) {
        if (inputW) {
          API::Run &run = inputW->mutableRun();
          // Check to see if LoadEventNexus had T0 from TOPAZ Parameter file
          IAlgorithm_sptr alg1 = createChildAlgorithm("ChangeBinOffset");
          alg1->setProperty<MatrixWorkspace_sptr>("InputWorkspace", inputW);
          alg1->setProperty<MatrixWorkspace_sptr>("OutputWorkspace", inputW);
          if (run.hasProperty("T0")) {
            double T0IDF = run.getPropertyValueAsType<double>("T0");
            alg1->setProperty("Offset", mT0 - T0IDF);
          } else {
            alg1->setProperty("Offset", mT0);
          }
          alg1->executeAsChildAlg();
          inputW = alg1->getProperty("OutputWorkspace");
          // set T0 in the run parameters
          run.addProperty<double>("T0", mT0, true);
        } else if (inputP) {
          // set T0 in the run parameters
          API::Run &run = inputP->mutableRun();
          run.addProperty<double>("T0", mT0, true);
234
        }
235
      }
236
    }
237

238
239
240
241
242
243
    if (line[0] != '5')
      continue;

    std::stringstream(line) >> count >> id >> nrows >> ncols >> width >>
        height >> depth >> detd >> x >> y >> z >> base_x >> base_y >> base_z >>
        up_x >> up_y >> up_z;
244
    if (id == 10 && filenames.size() == 2 && instname == "SNAP") {
245
      input.close();
246
      input.open(filenames[1].c_str());
247
248
249
250
251
252
253
254
255
      while (std::getline(input, line)) {
        if (line[0] != '5')
          continue;

        std::stringstream(line) >> count >> id >> nrows >> ncols >> width >>
            height >> depth >> detd >> x >> y >> z >> base_x >> base_y >>
            base_z >> up_x >> up_y >> up_z;
        if (id == 10)
          break;
256
      }
257
258
    }
    boost::shared_ptr<RectangularDetector> det;
259
    std::string bankName = getBankName(bankPart, id);
Hahn, Steven's avatar
Hahn, Steven committed
260
    auto matchingDetector = std::find_if(
Hahn, Steven's avatar
Hahn, Steven committed
261
        detList.begin(), detList.end(),
262
263
        [&bankName](const boost::shared_ptr<RectangularDetector> &detector) {
          return detector->getName() == bankName;
Hahn, Steven's avatar
Hahn, Steven committed
264
        });
Hahn, Steven's avatar
Hahn, Steven committed
265
266
    if (matchingDetector != detList.end()) {
      det = *matchingDetector;
Hahn, Steven's avatar
Hahn, Steven committed
267
    }
268
269
270
271

    V3D rX(base_x, base_y, base_z);
    V3D rY(up_x, up_y, up_z);

272
    if (det) {
273
      detname = det->getName();
274
      center(x, y, z, detname, ws, componentInfo);
275

276
277
278
      ComponentScaling detScaling;
      detScaling.scaleX = CM_TO_M * width / det->xsize();
      detScaling.scaleY = CM_TO_M * height / det->ysize();
279
      detScaling.componentName = detname;
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
      // Scaling will need both scale factors if LoadIsawPeaks or LoadIsawDetCal
      // has already
      // applied a calibration
      if (inputW) {
        Geometry::ParameterMap &pmap = inputW->instrumentParameters();
        auto oldscalex = pmap.getDouble(detname, std::string("scalex"));
        auto oldscaley = pmap.getDouble(detname, std::string("scaley"));
        if (!oldscalex.empty())
          detScaling.scaleX *= oldscalex[0];
        if (!oldscaley.empty())
          detScaling.scaleY *= oldscaley[0];
      }
      if (inputP) {
        Geometry::ParameterMap &pmap = inputP->instrumentParameters();
        auto oldscalex = pmap.getDouble(detname, std::string("scalex"));
        auto oldscaley = pmap.getDouble(detname, std::string("scaley"));
        if (!oldscalex.empty())
          detScaling.scaleX *= oldscalex[0];
        if (!oldscaley.empty())
          detScaling.scaleY *= oldscaley[0];
      }

302
      rectangularDetectorScalings.push_back(detScaling);
303

304
      doRotation(rX, rY, componentInfo, det);
305
    }
Hahn, Steven's avatar
Hahn, Steven committed
306
307
    auto bank = uniqueBanks.find(id);
    if (bank == uniqueBanks.end())
308
      continue;
Hahn, Steven's avatar
Hahn, Steven committed
309
    int idnum = *bank;
Hahn, Steven's avatar
Hahn, Steven committed
310

311
    bankName = getBankName(bankPart, idnum);
312
    // Retrieve it
Peterson, Peter's avatar
Peterson, Peter committed
313
    auto comp = inst->getComponentByName(bankName);
314
    // for Corelli with sixteenpack under bank
315
    if (instname == "CORELLI") {
316
317
      std::vector<Geometry::IComponent_const_sptr> children;
      boost::shared_ptr<const Geometry::ICompAssembly> asmb =
318
319
          boost::dynamic_pointer_cast<const Geometry::ICompAssembly>(
              inst->getComponentByName(bankName));
320
321
322
323
      asmb->getChildren(children, false);
      comp = children[0];
    }
    if (comp) {
324
      // Omitted scaling tubes
325
      detname = comp->getFullName();
326
      center(x, y, z, detname, ws, componentInfo);
327

328
329
      bool doWishCorrection =
          (instname == "WISH"); // TODO: find out why this is needed for WISH
330
      doRotation(rX, rY, componentInfo, comp, doWishCorrection);
331
    }
332
  }
333

334
  // Do this last, to avoid the issue of invalidating DetectorInfo
335
  applyScalings(ws, rectangularDetectorScalings);
336

337
  setProperty("InputWorkspace", ws);
338
}
339

340
341
342
343
344
345
346
347
/**
 * The intensity function calculates the intensity as a function of detector
 * position and angles
 * @param x :: The shift along the X-axis
 * @param y :: The shift along the Y-axis
 * @param z :: The shift along the Z-axis
 * @param detname :: The detector name
 * @param ws :: The workspace
348
 * @param componentInfo :: The component info object for the workspace
349
 */
Peterson, Peter's avatar
Peterson, Peter committed
350
void LoadIsawDetCal::center(const double x, const double y, const double z,
351
                            const std::string &detname, API::Workspace_sptr ws,
Owen Arnold's avatar
Owen Arnold committed
352
                            Geometry::ComponentInfo &componentInfo) {
353
354
355
356

  Instrument_sptr inst = getCheckInst(ws);

  IComponent_const_sptr comp = inst->getComponentByName(detname);
357
  if (comp == nullptr) {
358
359
    throw std::runtime_error("Component with name " + detname +
                             " was not found.");
360
361
  }

362
  const V3D position(x * CM_TO_M, y * CM_TO_M, z * CM_TO_M);
Peterson, Peter's avatar
Peterson, Peter committed
363

364
365
  const auto componentIndex = componentInfo.indexOf(comp->getComponentID());
  componentInfo.setPosition(componentIndex, position);
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
}

/**
 * Gets the instrument of the workspace, checking that the workspace
 * and the instrument are as expected.
 *
 * @param ws workspace (expected Matrix or Peaks Workspace)
 *
 * @throw std::runtime_error if there's any problem with the workspace or it is
 * not possible to get an instrument object from it
 */
Instrument_sptr LoadIsawDetCal::getCheckInst(API::Workspace_sptr ws) {
  MatrixWorkspace_sptr inputW =
      boost::dynamic_pointer_cast<MatrixWorkspace>(ws);
  PeaksWorkspace_sptr inputP = boost::dynamic_pointer_cast<PeaksWorkspace>(ws);

  // Get some stuff from the input workspace
  Instrument_sptr inst;
  if (inputW) {
    inst = boost::const_pointer_cast<Instrument>(inputW->getInstrument());
    if (!inst)
      throw std::runtime_error("Could not get a valid instrument from the "
                               "MatrixWorkspace provided as input");
  } else if (inputP) {
    inst = boost::const_pointer_cast<Instrument>(inputP->getInstrument());
    if (!inst)
      throw std::runtime_error("Could not get a valid instrument from the "
                               "PeaksWorkspace provided as input");
  } else {
Peterson, Peter's avatar
Peterson, Peter committed
395
396
397
398
    throw std::runtime_error("Could not get a valid instrument from the "
                             "workspace which does not seem to be valid as "
                             "input (must be either MatrixWorkspace or "
                             "PeaksWorkspace");
399
400
401
402
403
  }

  return inst;
}

404
405
406
407
std::vector<std::string> LoadIsawDetCal::getFilenames() {
  std::vector<std::string> filenamesFromPropertyUnraveld;
  std::vector<std::vector<std::string>> filenamesFromProperty =
      this->getProperty("Filename");
408
  for (const auto &outer : filenamesFromProperty) {
409
410
411
    std::copy(outer.begin(), outer.end(),
              std::back_inserter(filenamesFromPropertyUnraveld));
  }
412
413
414
415
416
417

  // shouldn't be used except for legacy cases
  const std::string filename2 = this->getProperty("Filename2");
  if (!filename2.empty())
    filenamesFromPropertyUnraveld.push_back(filename2);

418
419
420
  return filenamesFromPropertyUnraveld;
}

Ian Bush's avatar
Ian Bush committed
421
422
423
424
425
/**
 * Perform the rotation for the calibration
 *
 * @param rX the vector of (base_x, base_y, base_z) from the calibration file
 * @param rY the vector of (up_x, up_y, up_z) from the calibration file
426
 * @param componentInfo the ComponentInfo object from the workspace
Ian Bush's avatar
Ian Bush committed
427
428
429
 * @param comp the component to rotate
 * @param doWishCorrection if true apply a special correction for WISH
 */
430
void LoadIsawDetCal::doRotation(V3D rX, V3D rY, ComponentInfo &componentInfo,
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
                                boost::shared_ptr<const IComponent> comp,
                                bool doWishCorrection) {
  // These are the ISAW axes
  rX.normalize();
  rY.normalize();

  // These are the original axes
  const V3D oX(1., 0., 0.);
  const V3D oY(0., 1., 0.);

  // Axis that rotates X
  V3D ax1 = oX.cross_prod(rX);
  // Rotation angle from oX to rX
  double angle1 = oX.angle(rX) * DegreesPerRadian;
  if (doWishCorrection)
    angle1 += 180.0;
  // Create the first quaternion
  Quat Q1(angle1, ax1);

  // Now we rotate the original Y using Q1
  V3D roY = oY;
  Q1.rotate(roY);
  // Find the axis that rotates oYr onto rY
  V3D ax2 = roY.cross_prod(rY);
  const double angle2 = roY.angle(rY) * DegreesPerRadian;
  Quat Q2(angle2, ax2);

  // Final = those two rotations in succession; Q1 is done first.
  Quat Rot = Q2 * Q1;

  // Then find the corresponding relative position
462
463
464
  const auto componentIndex = componentInfo.indexOf(comp->getComponentID());

  componentInfo.setRotation(componentIndex, Rot);
465
466
}

Ian Bush's avatar
Ian Bush committed
467
468
469
470
471
/**
 * Apply the scalings from the calibration file. This is called after doing the
 *moves and rotations associated with the calibration, to avoid the problem of
 *invalidation DetectorInfo after writing to the parameter map.
 *
Ian Bush's avatar
Ian Bush committed
472
 * @param ws the input workspace
Ian Bush's avatar
Ian Bush committed
473
474
475
 * @param rectangularDetectorScalings a vector containing a component ID, and
 *values for scalex and scaley
 */
476
void LoadIsawDetCal::applyScalings(
477
    Workspace_sptr &ws,
478
479
480
    const std::vector<ComponentScaling> &rectangularDetectorScalings) {

  for (const auto &scaling : rectangularDetectorScalings) {
481
482
483
484
485
486
    IAlgorithm_sptr alg1 = createChildAlgorithm("ResizeRectangularDetector");
    alg1->setProperty<Workspace_sptr>("Workspace", ws);
    alg1->setProperty("ComponentName", scaling.componentName);
    alg1->setProperty("ScaleX", scaling.scaleX);
    alg1->setProperty("ScaleY", scaling.scaleY);
    alg1->executeAsChildAlg();
487
488
489
  }
}

490
491
} // namespace Algorithm
} // namespace Mantid