He3TubeEfficiency.cpp 18.2 KB
Newer Older
1
#include "MantidAlgorithms/He3TubeEfficiency.h"
2
3
4
#include "MantidAPI/HistogramValidator.h"
#include "MantidAPI/InstrumentValidator.h"
#include "MantidAPI/WorkspaceUnitValidator.h"
5
#include "MantidDataObjects/EventWorkspace.h"
6
#include "MantidDataObjects/Workspace2D.h"
7
#include "MantidKernel/ArrayBoundedValidator.h"
8
#include "MantidKernel/ArrayProperty.h"
9
10
#include "MantidKernel/CompositeValidator.h"

11
#include <algorithm>
12
#include <cmath>
13
#include <stdexcept>
14
15
16
17
18

// this should be a big number but not so big that there are rounding errors
const double DIST_TO_UNIVERSE_EDGE = 1e3;

// Scalar constant for exponential in units of K / (m Angstroms atm)
19
const double EXP_SCALAR_CONST = 2175.486863864;
20

21
22
23
// Tolerance for diameter/thickness comparison
const double TOL = 1.0e-8;

24
25
namespace Mantid {
namespace Algorithms {
26
27
28
29
// Register the class into the algorithm factory
DECLARE_ALGORITHM(He3TubeEfficiency)

/// Default constructor
30
31
32
He3TubeEfficiency::He3TubeEfficiency()
    : Algorithm(), inputWS(), outputWS(), paraMap(NULL), shapeCache(),
      samplePos(), spectraSkipped(), progress(NULL) {
33
  this->shapeCache.clear();
34
35
36
}

/// Destructor
37
38
He3TubeEfficiency::~He3TubeEfficiency() {
  if (this->progress) {
39
40
    delete this->progress;
  }
41
42
43
44
45
}

/**
 * Declare algorithm properties
 */
46
void He3TubeEfficiency::init() {
47
48
49
50
51
52
  using namespace Mantid::Kernel;

  auto wsValidator = boost::make_shared<CompositeValidator>();
  wsValidator->add<API::WorkspaceUnitValidator>("Wavelength");
  wsValidator->add<API::HistogramValidator>();
  wsValidator->add<API::InstrumentValidator>();
53
54
55
56
57
58
59
60
61
  this->declareProperty(
      new API::WorkspaceProperty<API::MatrixWorkspace>(
          "InputWorkspace", "", Kernel::Direction::Input, wsValidator),
      "Name of the input workspace");
  this->declareProperty(
      new API::WorkspaceProperty<API::MatrixWorkspace>(
          "OutputWorkspace", "", Kernel::Direction::Output),
      "Name of the output workspace, can be the same as the input");
  auto mustBePositive = boost::make_shared<Kernel::BoundedValidator<double>>();
62
  mustBePositive->setLower(0.0);
63
64
65
  this->declareProperty(
      new Kernel::PropertyWithValue<double>("ScaleFactor", 1.0, mustBePositive),
      "Constant factor with which to scale the calculated"
66
      "detector efficiency. Same factor applies to all efficiencies.");
67

68
69
  auto mustBePosArr =
      boost::make_shared<Kernel::ArrayBoundedValidator<double>>();
70
  mustBePosArr->setLower(0.0);
71
72
  this->declareProperty(
      new Kernel::ArrayProperty<double>("TubePressure", mustBePosArr),
73
74
      "Provide overriding the default tube pressure. The pressure must "
      "be specified in atm.");
75
76
  this->declareProperty(
      new Kernel::ArrayProperty<double>("TubeThickness", mustBePosArr),
77
78
      "Provide overriding the default tube thickness. The thickness must "
      "be specified in metres.");
79
80
  this->declareProperty(
      new Kernel::ArrayProperty<double>("TubeTemperature", mustBePosArr),
81
82
      "Provide overriding the default tube temperature. The temperature must "
      "be specified in Kelvin.");
83
84
85
86
87
}

/**
 * Executes the algorithm
 */
88
void He3TubeEfficiency::exec() {
89
90
  // Get the workspaces
  this->inputWS = this->getProperty("InputWorkspace");
91
  this->outputWS = this->getProperty("OutputWorkspace");
92

93
  if (this->outputWS != this->inputWS) {
94
    this->outputWS = API::WorkspaceFactory::Instance().create(this->inputWS);
95
  }
96

97
98
99
  // Get the detector parameters
  this->paraMap = &(this->inputWS->instrumentParameters());

100
101
  // Store some information about the instrument setup that will not change
  this->samplePos = this->inputWS->getInstrument()->getSample()->getPos();
102

103
  // Check if it is an event workspace
104
105
106
  DataObjects::EventWorkspace_const_sptr eventW =
      boost::dynamic_pointer_cast<const DataObjects::EventWorkspace>(inputWS);
  if (eventW != NULL) {
107
108
    this->execEvent();
    return;
109
110
  }

111
  std::size_t numHists = this->inputWS->getNumberHistograms();
112
  this->progress = new API::Progress(this, 0.0, 1.0, numHists);
113
114

  PARALLEL_FOR2(inputWS, outputWS)
115
  for (int i = 0; i < static_cast<int>(numHists); ++i) {
116
117
118
    PARALLEL_START_INTERUPT_REGION

    this->outputWS->setX(i, this->inputWS->refX(i));
119
    try {
120
      this->correctForEfficiency(i);
121
    } catch (std::out_of_range &) {
122
123
      // if we don't have all the data there will be spectra we can't correct,
      // avoid leaving the workspace part corrected
124
      Mantid::MantidVec &dud = this->outputWS->dataY(i);
125
      std::transform(dud.begin(), dud.end(), dud.begin(),
126
127
                     std::bind2nd(std::multiplies<double>(), 0));
      PARALLEL_CRITICAL(deteff_invalid) {
128
129
130
131
        this->spectraSkipped.push_back(this->inputWS->getAxis(1)->spectraNo(i));
      }
    }

132
133
134
    // make regular progress reports
    this->progress->report();
    // check for canceling the algorithm
135
    if (i % 1000 == 0) {
136
137
138
139
140
141
142
143
144
145
146
147
      interruption_point();
    }

    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  this->logErrors();
  this->setProperty("OutputWorkspace", this->outputWS);
}

/**
148
149
150
 * Corrects a spectra for the detector efficiency calculated from detector
 * information. Gets the detector information and uses this to calculate its
 * efficiency
151
 *  @param spectraIndex :: index of the spectrum to get the efficiency for
152
153
154
155
156
 *  @throw invalid_argument if the shape of a detector is isn't a cylinder
 *  aligned along one axis
 *  @throw NotFoundError if the detector or its gas pressure or wall thickness
 *  were not found
 */
157
void He3TubeEfficiency::correctForEfficiency(std::size_t spectraIndex) {
158
  Geometry::IDetector_const_sptr det = this->inputWS->getDetector(spectraIndex);
159
  if (det->isMonitor() || det->isMasked()) {
160
161
162
    return;
  }

163
164
165
  const double exp_constant = this->calculateExponential(spectraIndex, det);
  const double scale = this->getProperty("ScaleFactor");

166
167
168
169
170
171
  Mantid::MantidVec &yout = this->outputWS->dataY(spectraIndex);
  Mantid::MantidVec &eout = this->outputWS->dataE(spectraIndex);
  // Need the original values so this is not a reference
  const Mantid::MantidVec yValues = this->inputWS->readY(spectraIndex);
  const Mantid::MantidVec eValues = this->inputWS->readE(spectraIndex);

172
173
174
  auto yinItr = yValues.cbegin();
  auto einItr = eValues.cbegin();
  auto xItr = this->inputWS->readX(spectraIndex).cbegin();
175
176
  auto youtItr = yout.begin();
  auto eoutItr = eout.begin();
177

178
  for (; youtItr != yout.end(); ++youtItr, ++eoutItr) {
179
    const double wavelength = (*xItr + *(xItr + 1)) / 2.0;
180
181
    const double effcorr =
        this->detectorEfficiency(exp_constant * wavelength, scale);
182
183
    *youtItr = (*yinItr) * effcorr;
    *eoutItr = (*einItr) * effcorr;
184
185
    ++yinItr;
    ++einItr;
186
187
188
189
190
191
192
193
194
195
196
    ++xItr;
  }

  return;
}

/**
 * This function calculates the exponential contribution to the He3 tube
 * efficiency.
 * @param spectraIndex :: the current index to calculate
 * @param idet :: the current detector pointer
197
 * @throw out_of_range if twice tube thickness is greater than tube diameter
198
199
 * @return the exponential contribution for the given detector
 */
200
201
202
double He3TubeEfficiency::calculateExponential(
    std::size_t spectraIndex,
    boost::shared_ptr<const Geometry::IDetector> idet) {
203
  // Get the parameters for the current associated tube
204
205
206
207
  double pressure =
      this->getParameter("TubePressure", spectraIndex, "tube_pressure", idet);
  double tubethickness =
      this->getParameter("TubeThickness", spectraIndex, "tube_thickness", idet);
208
  double temperature = this->getParameter("TubeTemperature", spectraIndex,
209
                                          "tube_temperature", idet);
210
211

  double detRadius(0.0);
212
  Kernel::V3D detAxis;
213
  this->getDetectorGeometry(idet, detRadius, detAxis);
214
215
  double detDiameter = 2.0 * detRadius;
  double twiceTubeThickness = 2.0 * tubethickness;
216
217
218
219

  // now get the sin of the angle, it's the magnitude of the cross product of
  // unit vector along the detector tube axis and a unit vector directed from
  // the sample to the detector center
220
  Kernel::V3D vectorFromSample = idet->getPos() - this->samplePos;
221
  vectorFromSample.normalize();
222
  Kernel::Quat rot = idet->getRotation();
223
224
225
226
227
228
  // rotate the original cylinder object axis to get the detector axis in the
  // actual instrument
  rot.rotate(detAxis);
  detAxis.normalize();
  // Scalar product is quicker than cross product
  double cosTheta = detAxis.scalar_prod(vectorFromSample);
229
  double sinTheta = std::sqrt(1.0 - cosTheta * cosTheta);
230

231
  const double straight_path = detDiameter - twiceTubeThickness;
232
233
234
  if (std::fabs(straight_path - 0.0) < TOL) {
    throw std::out_of_range("Twice tube thickness cannot be greater than "
                            "or equal to the tube diameter");
235
236
237
  }

  const double pathlength = straight_path / sinTheta;
238
  return EXP_SCALAR_CONST * (pressure / temperature) * pathlength;
239
240
}

241
242
/**
 * Update the shape cache if necessary
243
244
245
 * @param det :: a pointer to the detector to query
 * @param detRadius :: An output parameter that contains the detector radius
 * @param detAxis :: An output parameter that contains the detector axis vector
246
 */
247
void He3TubeEfficiency::getDetectorGeometry(
248
    const boost::shared_ptr<const Geometry::IDetector> &det, double &detRadius,
249
    Kernel::V3D &detAxis) {
250
  boost::shared_ptr<const Geometry::Object> shape_sptr = det->shape();
251
252
253
254
255
256
257
  if (!shape_sptr) {
    throw std::runtime_error(
        "Detector geometry error: detector with id: " +
        boost::lexical_cast<std::string>(det->getID()) +
        " does not have shape. Is this a detectors group?\n"
        "The algorithm works for instruments with one-to-one "
        "spectra-to-detector maps only!");
258
  }
259
260
261
262
263
264
265
266
267
  std::map<const Geometry::Object *,
           std::pair<double, Kernel::V3D>>::const_iterator it =
      this->shapeCache.find(shape_sptr.get());
  if (it == this->shapeCache.end()) {
    double xDist = distToSurface(Kernel::V3D(DIST_TO_UNIVERSE_EDGE, 0, 0),
                                 shape_sptr.get());
    double zDist = distToSurface(Kernel::V3D(0, 0, DIST_TO_UNIVERSE_EDGE),
                                 shape_sptr.get());
    if (std::abs(zDist - xDist) < 1e-8) {
268
      detRadius = zDist / 2.0;
269
      detAxis = Kernel::V3D(0, 1, 0);
270
      // assume radii in z and x and the axis is in the y
271
272
273
274
      PARALLEL_CRITICAL(deteff_shapecachea) {
        this->shapeCache.insert(
            std::pair<const Geometry::Object *, std::pair<double, Kernel::V3D>>(
                shape_sptr.get(),
275
                std::pair<double, Kernel::V3D>(detRadius, detAxis)));
276
277
278
      }
      return;
    }
279
280
281
    double yDist = distToSurface(Kernel::V3D(0, DIST_TO_UNIVERSE_EDGE, 0),
                                 shape_sptr.get());
    if (std::abs(yDist - zDist) < 1e-8) {
282
      detRadius = yDist / 2.0;
283
      detAxis = Kernel::V3D(1, 0, 0);
284
      // assume that y and z are radii of the cylinder's circular cross-section
285
      // and the axis is perpendicular, in the x direction
286
287
288
289
      PARALLEL_CRITICAL(deteff_shapecacheb) {
        this->shapeCache.insert(
            std::pair<const Geometry::Object *, std::pair<double, Kernel::V3D>>(
                shape_sptr.get(),
290
                std::pair<double, Kernel::V3D>(detRadius, detAxis)));
291
292
293
294
      }
      return;
    }

295
    if (std::abs(xDist - yDist) < 1e-8) {
296
      detRadius = xDist / 2.0;
297
      detAxis = Kernel::V3D(0, 0, 1);
298
299
300
301
      PARALLEL_CRITICAL(deteff_shapecachec) {
        this->shapeCache.insert(
            std::pair<const Geometry::Object *, std::pair<double, Kernel::V3D>>(
                shape_sptr.get(),
302
                std::pair<double, Kernel::V3D>(detRadius, detAxis)));
303
304
305
      }
      return;
    }
306
  } else {
307
    std::pair<double, Kernel::V3D> geometry = it->second;
308
309
    detRadius = geometry.first;
    detAxis = geometry.second;
310
311
312
313
314
315
  }
}

/**
 * For basic shapes centered on the origin (0,0,0) this returns the distance to
 * the surface in the direction of the point given
316
 *  @param start :: the distance calculated from origin to the surface in a line
317
 *  towards this point. It should be outside the shape
318
319
 *  @param shape :: the object to calculate for, should be centered on the
 * origin
320
321
322
323
 *  @return the distance to the surface in the direction of the point given
 *  @throw invalid_argument if there is any error finding the distance
 * @returns The distance to the surface in metres
 */
324
double He3TubeEfficiency::distToSurface(const Kernel::V3D start,
325
                                        const Geometry::Object *shape) const {
326
  // get a vector from the point that was passed to the origin
327
  Kernel::V3D direction = Kernel::V3D(0.0, 0.0, 0.0) - start;
328
329
330
331
  // it needs to be a unit vector
  direction.normalize();
  // put the point and the vector (direction) together to get a line,
  // here called a track
332
  Geometry::Track track(start, direction);
333
334
335
336
  // split the track (line) up into the part that is inside the shape and the
  // part that is outside
  shape->interceptSurface(track);

337
  if (track.count() != 1) {
338
339
    // the track missed the shape, probably the shape is not centered on
    // the origin
340
341
    throw std::invalid_argument(
        "Fatal error interpreting the shape of a detector");
342
343
344
  }
  // the first part of the track will be the part inside the shape,
  // return its length
345
  return track.cbegin()->distInsideObject;
346
347
}

348
349
350
/**
 * Calculate the detector efficiency from the detector parameters and the
 * spectrum's x-axis.
351
352
 * @param alpha :: the value to feed to the exponential
 * @param scale_factor :: an overall value for scaling the efficiency
353
 * @return the calculated efficiency
354
 */
355
double He3TubeEfficiency::detectorEfficiency(const double alpha,
356
                                             const double scale_factor) const {
357
  return (scale_factor / (1.0 - std::exp(-alpha)));
358
359
360
361
362
}

/**
 * Logs if there were any problems locating spectra.
 */
363
void He3TubeEfficiency::logErrors() const {
364
  std::vector<int>::size_type nspecs = this->spectraSkipped.size();
365
366
367
  if (nspecs > 0) {
    this->g_log.warning() << "There were " << nspecs
                          << " spectra that could not be corrected. ";
368
    this->g_log.debug() << "Unaffected spectra numbers: ";
369
    for (size_t i = 0; i < nspecs; ++i) {
370
      this->g_log.debug() << this->spectraSkipped[i] << " ";
371
    }
372
    this->g_log.debug() << std::endl;
373
374
  }
}
375
376
377
378

/**
 * Retrieve the detector parameter either from the workspace property or from
 * the associated detector property.
379
380
381
382
 * @param wsPropName :: the workspace property name for the detector parameter
 * @param currentIndex :: the currently requested spectra index
 * @param detPropName :: the detector property name for the detector parameter
 * @param idet :: the current detector
383
384
 * @return the value of the detector property
 */
385
386
387
double He3TubeEfficiency::getParameter(
    std::string wsPropName, std::size_t currentIndex, std::string detPropName,
    boost::shared_ptr<const Geometry::IDetector> idet) {
388
389
  std::vector<double> wsProp = this->getProperty(wsPropName);

390
  if (wsProp.empty()) {
391
    return idet->getNumberParameter(detPropName).at(0);
392
393
  } else {
    if (wsProp.size() == 1) {
394
      return wsProp.at(0);
395
    } else {
396
397
398
399
      return wsProp.at(currentIndex);
    }
  }
}
400
401
402
403

/**
 * Execute for events
 */
404
void He3TubeEfficiency::execEvent() {
405
  this->g_log.information("Processing event workspace");
406

407
408
409
410
411
  const API::MatrixWorkspace_const_sptr matrixInputWS =
      this->getProperty("InputWorkspace");
  DataObjects::EventWorkspace_const_sptr inputWS =
      boost::dynamic_pointer_cast<const DataObjects::EventWorkspace>(
          matrixInputWS);
412
413

  // generate the output workspace pointer
414
415
  API::MatrixWorkspace_sptr matrixOutputWS =
      this->getProperty("OutputWorkspace");
416
  DataObjects::EventWorkspace_sptr outputWS;
417
418
419
420
  if (matrixOutputWS == matrixInputWS) {
    outputWS = boost::dynamic_pointer_cast<DataObjects::EventWorkspace>(
        matrixOutputWS);
  } else {
421
422
    // Make a brand new EventWorkspace
    outputWS = boost::dynamic_pointer_cast<DataObjects::EventWorkspace>(
423
424
        API::WorkspaceFactory::Instance().create(
            "EventWorkspace", inputWS->getNumberHistograms(), 2, 1));
425
    // Copy geometry over.
426
427
    API::WorkspaceFactory::Instance().initializeFromParent(inputWS, outputWS,
                                                           false);
428
    // You need to copy over the data as well.
429
    outputWS->copyDataFrom((*inputWS));
430
431

    // Cast to the matrixOutputWS and save it
432
433
    matrixOutputWS =
        boost::dynamic_pointer_cast<API::MatrixWorkspace>(outputWS);
434
435
436
    this->setProperty("OutputWorkspace", matrixOutputWS);
  }

437
  std::size_t numHistograms = inputWS->getNumberHistograms();
438
439
  this->progress = new API::Progress(this, 0.0, 1.0, numHistograms);
  PARALLEL_FOR1(outputWS)
440
  for (int i = 0; i < static_cast<int>(numHistograms); ++i) {
441
442
    PARALLEL_START_INTERUPT_REGION

443
    Geometry::IDetector_const_sptr det = inputWS->getDetector(i);
444
    if (det->isMonitor() || det->isMasked()) {
445
446
447
      continue;
    }

448
    double exp_constant = 0.0;
449
    try {
450
      exp_constant = this->calculateExponential(i, det);
451
    } catch (std::out_of_range &) {
452
      // Parameters are bad so skip correction
453
      PARALLEL_CRITICAL(deteff_invalid) {
454
        this->spectraSkipped.push_back(inputWS->getAxis(1)->spectraNo(i));
455
        outputWS->maskWorkspaceIndex(i);
456
      }
457
458
    }

459
460
    // Do the correction
    DataObjects::EventList *evlist = outputWS->getEventListPtr(i);
461
    switch (evlist->getEventType()) {
462
463
464
    case API::TOF:
      // Switch to weights if needed.
      evlist->switchTo(API::WEIGHTED);
465
    // Fall through
466
    case API::WEIGHTED:
467
      eventHelper(evlist->getWeightedEvents(), exp_constant);
468
469
      break;
    case API::WEIGHTED_NOTIME:
470
      eventHelper(evlist->getWeightedEventsNoTime(), exp_constant);
471
472
473
474
475
476
      break;
    }

    this->progress->report();

    // check for canceling the algorithm
477
    if (i % 1000 == 0) {
478
479
480
481
482
483
484
485
      interruption_point();
    }

    PARALLEL_END_INTERUPT_REGION
  }
  PARALLEL_CHECK_INTERUPT_REGION

  outputWS->clearMRU();
486
487

  this->logErrors();
488
489
}

490
491
492
/**
 * Private function for doing the event correction.
 * @param events :: the list of events to correct
493
 * @param expval :: the value of the exponent for the detector efficiency
494
 */
495
496
template <class T>
void He3TubeEfficiency::eventHelper(std::vector<T> &events, double expval) {
497
  const double scale = this->getProperty("ScaleFactor");
498
  typename std::vector<T>::iterator it;
499
500
501
  for (it = events.begin(); it != events.end(); ++it) {
    float de =
        static_cast<float>(this->detectorEfficiency(expval * it->tof(), scale));
502
    it->m_weight *= de;
503
    it->m_errorSquared *= de * de;
504
505
506
  }
}

507
508
} // namespace Algorithms
} // namespace Mantid