68using namespace Kernel;
69using namespace Geometry;
71using namespace DataObjects;
75 :
API::
Algorithm(), m_smu(0.), m_amu(0.), m_radius(0.), m_power_th(0.), m_lamda_weight(),
76 m_onlySphericalAbsorption(false), m_returnTransmissionOnly(false), m_useScaleFactors(false) {}
81 auto wsValidator = std::make_shared<InstrumentValidator>();
85 "The X values for the input workspace must be in units of wavelength or TOF");
87 "Output workspace name");
89 auto mustBePositive = std::make_shared<BoundedValidator<double>>();
90 mustBePositive->setLower(0.0);
92 "Linear scattering coefficient in 1/cm. "
93 "If not provided this will be calculated from the "
94 "material cross-section if present (set with SetSampleMaterial)");
96 "Linear absorption coefficient at 1.8 Angstroms in 1/cm. "
97 "If not provided this will be calculated from the "
98 "material cross-section if present (set with SetSampleMaterial)");
100 "Radius of the sample in centimeters. f not provided the "
101 "radius will be taken from the sample shape if it is a sphere "
102 "(set with SetSample).");
104 "Keep the output workspace as an EventWorkspace, if the "
105 "input has events (default).\n"
106 "If false, then the workspace gets converted to a "
107 "Workspace2D histogram.");
109 "All corrections done if false (default).\n"
110 "If true, only the spherical absorption correction.");
112 "Corrections applied to data if false (default).\n"
113 "If true, only return the transmission coefficient.");
116 "No scale factors if false (default).\n"
117 "If true, use scale factors from instrument parameter map.");
123 std::map<std::string, std::string> result;
129 if (!inputWorkspace) {
130 result[
"InputWorkspace"] =
"The InputWorkspace must be a MatrixWorkspace.";
134 const auto &sampleShape = inputWorkspace->sample().getShape();
135 if (!sampleShape.hasValidShape() ||
137 result[
"Radius"] =
"Please supply a radius or provide a workspace with a spherical sample set.";
155 g_log.
warning() <<
"Lorentz Correction was already done for this "
156 "workspace. OnlySphericalAbsorption was changed to "
162 const std::string &unitStr =
m_inputWS->getAxis(0)->unit()->unitID();
183 const auto numHists =
static_cast<int64_t
>(
m_inputWS->getNumberHistograms());
184 const auto specSize =
static_cast<int64_t
>(
m_inputWS->blocksize());
186 throw std::runtime_error(
"Problem in AnvredCorrection::events not binned");
189 const auto &spectrumInfo =
m_inputWS->spectrumInfo();
190 double L1 = spectrumInfo.l1();
192 Progress prog(
this, 0.0, 1.0, numHists);
195 for (int64_t i = 0; i < int64_t(numHists); ++i) {
199 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
207 double L2 =
pmap.at(UnitParams::l2);
208 double scattering =
pmap.at(UnitParams::twoTheta);
212 double pathlength = 0.0;
214 std::string bankName;
216 const auto &det = spectrumInfo.detector(i);
217 bankName = det.getParent()->getParent()->getName();
218 scale_init(inst, L2, depth, pathlength, bankName);
224 const auto &inSpec =
m_inputWS->getSpectrum(i);
225 correctionFactors->setSharedX(i, inSpec.sharedX());
228 const auto &Yin = inSpec.y();
229 const auto &Ein = inSpec.x();
232 auto &
Y = correctionFactors->mutableY(i);
233 auto &E = correctionFactors->mutableE(i);
235 bool muRTooLarge =
false;
236 for (int64_t j = 0; j < specSize; j++) {
243 const auto eventWeight = this->
getEventWeight(lambda, scattering, muRTooLarge);
244 double scaleFactor(eventWeight);
246 scaleFactor =
scale_exec(bankName,
lambda, depth, inst, pathlength, eventWeight);
248 Y[j] = Yin[j] * scaleFactor;
249 E[j] = Ein[j] * scaleFactor;
254 g_log.
warning(
"Absorption correction not accurate for muR > 8 which was exceeded in spectrum index " +
265 API::Run &run = correctionFactors->mutableRun();
268 run.addProperty<
bool>(
"LorentzCorrection",
true,
true);
279 const auto numHists =
static_cast<int64_t
>(
m_inputWS->getNumberHistograms());
280 std::string unitStr =
m_inputWS->getAxis(0)->unit()->unitID();
281 auto correctionFactors = create<EventWorkspace>(*
m_inputWS);
282 correctionFactors->sortAll(
TOF_SORT,
nullptr);
285 g_log.
debug(
"Correcting EventWorkspace in-place.");
290 const auto &spectrumInfo =
eventW->spectrumInfo();
291 double L1 = spectrumInfo.l1();
293 Progress prog(
this, 0.0, 1.0, numHists);
296 for (int64_t i = 0; i < int64_t(numHists); ++i) {
300 correctionFactors->setHistogram(i,
eventW->binEdges(i));
303 if (!spectrumInfo.hasDetectors(i) || spectrumInfo.isMonitor(i))
310 double L2 =
pmap.at(UnitParams::l2);
311 double scattering =
pmap.at(UnitParams::twoTheta);
318 double pathlength = 0.0;
319 std::string bankName;
321 const auto &det = spectrumInfo.detector(i);
322 bankName = det.getParent()->getParent()->getName();
323 scale_init(inst, L2, depth, pathlength, bankName);
327 bool muRTooLarge =
false;
329 for (
auto &ev : events) {
333 if (
"TOF" == unitStr)
341 ev.m_errorSquared =
static_cast<float>(ev.m_errorSquared *
value *
value);
342 ev.m_weight *=
static_cast<float>(
value);
346 g_log.
warning(
"Absorption correction not accurate for muR > 9 cm^-1 which was exceeded in spectrum index " +
350 correctionFactors->getSpectrum(i) += events;
354 eventW->getSpectrum(i).clear();
363 API::Run &run = correctionFactors->mutableRun();
366 run.addProperty<
bool>(
"LorentzCorrection",
true,
true);
367 setProperty(
"OutputWorkspace", std::move(correctionFactors));
383 if (scatterXSection != 0.0) {
391 m_inputWS->sample().getShape().shapeInfo().sphereGeometry().radius * 100;
395 auto shape = std::shared_ptr<IObject>(
396 m_inputWS->sample().getShape().cloneWithMaterial(
Material(
"SetInAnvredCorrection", neutron, 1.0)));
397 m_inputWS->mutableSample().setShape(shape);
401 <<
"LinearAbsorptionCoef = " <<
m_amu <<
" 1/cm\n"
402 <<
"Radius = " <<
m_radius <<
" cm\n"
403 <<
"Power Lorentz corrections = " <<
m_power_th <<
" \n";
433 double sin_theta = std::sin(two_theta / 2);
434 double pix_weight = sin_theta * sin_theta;
436 double event_weight = pix_weight * lamda_w * transinv;
470 throw std::runtime_error(
"muR cannot be negative");
471 }
else if (mur > 8.0) {
475 auto theta = 0.5 * twoth *
radtodeg;
476 if (theta < 0. || theta > 90.) {
477 std::ostringstream s;
479 throw std::runtime_error(
"theta is not in allowed range :" + s.str());
495 auto ith =
static_cast<size_t>(theta / 5.);
498 size_t ncoef =
sizeof pc /
sizeof pc[0];
499 for (
size_t icoef = 0; icoef < ncoef; icoef++) {
500 lnA_1 = lnA_1 * mur +
pc[icoef][ith];
501 lnA_2 = lnA_2 * mur +
pc[icoef][ith + 1];
503 double A1 = std::exp(lnA_1);
504 double sin_th1_sq = std::pow(sin((
static_cast<double>(ith) * 5.0) /
radtodeg), 2);
505 double A2 = std::exp(lnA_2);
506 double sin_th2_sq = std::pow(sin((
static_cast<double>(ith + 1) * 5.0) /
radtodeg), 2);
509 double L1 = (A1 - A2) / (sin_th1_sq - sin_th2_sq);
510 double L0 = A1 - L1 * sin_th1_sq;
513 return 1 / (L0 + L1 * std::pow(sin(theta /
radtodeg), 2));
551 double &pathlength,
const std::string &bankName) {
553 std::shared_ptr<const IComponent> det0 = inst->getComponentByName(bankName);
554 if (
"CORELLI" == inst->getName())
556 std::vector<Geometry::IComponent_const_sptr> children;
557 auto asmb = std::dynamic_pointer_cast<const Geometry::ICompAssembly>(inst->getComponentByName(bankName));
558 asmb->getChildren(children,
false);
562 double cosA = det0->getDistance(*sample) / L2;
563 pathlength = depth / cosA;
569 const double mu = (9.614 *
lambda) + 0.266;
570 const double eff_center = 1.0 - std::exp(-
mu * depth);
571 const double eff_R = 1.0 - exp(-
mu * pathlength);
572 double scaleFactor(eventWeight * eff_center / eff_R);
574 bankName.erase(remove_if(bankName.begin(), bankName.end(), std::not_fn(::isdigit)), bankName.end());
575 if (inst->hasParameter(
"detScale" + bankName))
576 scaleFactor *=
static_cast<double>(inst->getNumberParameter(
"detScale" + bankName)[0]);
#define DECLARE_ALGORITHM(classname)
const std::vector< double > * lambda
double value
The value of the point.
#define PARALLEL_START_INTERRUPT_REGION
Begins a block to skip processing is the algorithm has been interupted Note the end of the block if n...
#define PARALLEL_END_INTERRUPT_REGION
Ends a block to skip processing is the algorithm has been interupted Note the start of the block if n...
#define PARALLEL_FOR_IF(condition)
Empty definitions - to enable set your complier to enable openMP.
#define PARALLEL_CHECK_INTERRUPT_REGION
Adds a check after a Parallel region to see if it was interupted.
Base class from which all concrete algorithm classes should be derived.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
std::string getPropertyValue(const std::string &name) const override
Get the value of a property as a string.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
bool hasProperty(const std::string &name) const
Does the property exist on the object.
void addProperty(Kernel::Property *prop, bool overwrite=false)
Add data to the object in the form of a property.
HeldType getPropertyValueAsType(const std::string &name) const
Get the value of a property as the given TYPE.
Helper class for reporting progress from algorithms.
This class stores information regarding an experimental run as a series of log entries.
A property class for workspaces.
double m_smu
linear scattering coefficient in 1/cm
virtual void defineProperties()
A virtual function in which additional properties of an algorithm should be declared.
double m_radius
sample radius in cm
bool m_onlySphericalAbsorption
std::map< std::string, std::string > validateInputs() override
validate inputs
void execEvent()
Event execution code.
virtual void retrieveProperties()
A virtual function in which additional properties should be retrieved into member variables.
double getEventWeight(const double lamda, const double two_theta, bool &muRTooLarge)
Get the weight factor that would be used for an event occuring at the specified wavelength,...
void scale_init(const Geometry::Instrument_const_sptr &inst, const double L2, const double depth, double &pathlength, const std::string &bankName)
bool m_returnTransmissionOnly
void init() override
Initialisation code.
double m_amu
linear absoprtion coefficient in 1/cm
void cleanup()
Algorithm cleanup.
void exec() override
Execution code.
double absor_sphere(const double twoth, const double wl, bool &muRTooLarge)
function to calculate a spherical absorption correction and tbar.
DataObjects::EventWorkspace_sptr eventW
Shared pointer to the event workspace.
std::vector< double > m_lamda_weight
lmabda weights
double m_power_th
Power of lamda in BuildLamdaWeights.
static double calc_Astar(const double theta, const double mur)
void BuildLamdaWeights()
Build the list of weights corresponding to different wavelengths.
AnvredCorrection()
(Empty) Constructor
API::MatrixWorkspace_sptr m_inputWS
A pointer to the input workspace.
double scale_exec(std::string &bankName, const double lambda, const double depth, const Geometry::Instrument_const_sptr &inst, const double pathlength, double eventWeight)
void retrieveBaseProperties()
Fetch the properties and set the appropriate member variables.
std::vector< WeightedEventNoTime > & getWeightedEventsNoTime()
Return the list of WeightedEvent contained.
void switchTo(Mantid::API::EventType newType) override
Switch the EventList to use the given EventType (TOF, WEIGHTED, or WEIGHTED_NOTIME)
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void debug(const std::string &msg)
Logs at debug level.
void notice(const std::string &msg)
Logs at notice level.
void warning(const std::string &msg)
Logs at warning level.
A material is defined as being composed of a given element, defined as a PhysicalConstants::NeutronAt...
double numberDensity() const
Get the number density.
double absorbXSection(const double lambda=PhysicalConstants::NeutronAtom::ReferenceLambda) const
Get the absorption cross section at a given wavelength in barns.
double totalScatterXSection() const
Return the total scattering cross section for a given wavelength in barns.
void report()
Increments the loop counter by 1, then sends the progress notification on behalf of its algorithm.
static T & Instance()
Return a reference to the Singleton instance, creating it if it does not already exist Creation is do...
double convertSingleFromTOF(const double xvalue, const double &l1, const int &emode, const UnitParametersMap ¶ms)
Convert from the time-of-flight to the concrete unit.
Time of flight in microseconds.
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
const int NUM_WAVELENGTHS
const double STEPS_PER_ANGSTROM
std::shared_ptr< const IComponent > IComponent_const_sptr
Typdef of a shared pointer to a const IComponent.
std::shared_ptr< const Instrument > Instrument_const_sptr
Shared pointer to an const instrument object.
std::unordered_map< UnitParams, double > UnitParametersMap
std::enable_if< std::is_pointer< Arg >::value, bool >::type threadSafe(Arg workspace)
Thread-safety check Checks the workspace to ensure it is suitable for multithreaded access.
A namespace containing physical constants that are required by algorithms and unit routines.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Generate a tableworkspace to store the calibration results.
std::string to_string(const wide_integer< Bits, Signed > &n)
@ Input
An input workspace.
@ Output
An output workspace.
Structure to store neutronic scattering information for the various elements.
static const double ReferenceLambda
The reference wavelength value for absorption cross sections.