42std::pair<double, double> calculateSummationTerm(
const Kernel::Material &material) {
43 std::pair<double, double> s;
45 const auto &formula = material.chemicalFormula();
46 auto sumLambda_first = [](
double sum,
const auto &formula_unit) {
47 return sum + formula_unit.multiplicity * formula_unit.atom->neutron.tot_scatt_xs / formula_unit.atom->mass;
49 auto sumLambda_second = [](
double sum,
const auto &formula_unit) {
50 const double mass_sq = formula_unit.atom->mass * formula_unit.atom->mass;
51 return sum + formula_unit.multiplicity * formula_unit.atom->neutron.tot_scatt_xs / mass_sq;
53 const double unnormalizedTermFirst = std::accumulate(formula.begin(), formula.end(), 0.0, sumLambda_first);
54 const double unnormalizedTermSecond = std::accumulate(formula.begin(), formula.end(), 0.0, sumLambda_second);
58 constexpr double neutronMassSq = neutronMass * neutronMass;
60 const double totalStoich = material.totalAtoms();
62 s.first = neutronMass * unnormalizedTermFirst / (4. * M_PI * totalStoich);
63 s.second = neutronMassSq * unnormalizedTermSecond / (4. * M_PI * totalStoich);
87 return "Calculate 1st or 2nd order Placzek correction factors using given workspace and incident spectrums.";
92 return {
"CalculatePlaczekSelfScattering",
"He3TubeEfficiency"};
103 auto wsValidator = std::make_shared<Mantid::Kernel::CompositeValidator>();
108 "Raw diffraction data workspace for associated correction to be "
109 "calculated for. Workspace must have instrument and sample data.");
111 auto inspValidator = std::make_shared<Mantid::Kernel::CompositeValidator>();
115 "Workspace of fitted incident spectrum with its derivatives (1st &| 2nd).");
120 "Workspace of efficiency spectrum with its derivatives (1st &| 2nd)."
121 "Default (not specified) will use LambdaD to calculate the efficiency spectrum.");
122 auto lambdadValidator = std::make_shared<Kernel::BoundedValidator<double>>();
123 lambdadValidator->setExclusive(
true);
124 lambdadValidator->setLower(0.0);
126 "LambdaD", 1.44, lambdadValidator,
127 "Reference wavelength in Angstrom, related to detector efficient coefficient alpha."
128 "The coefficient used to generate a generic detector efficiency curve,"
129 "eps = 1 - exp(1 - alpha*lambda), where alpha is 1/LambdaD."
130 "Default is set to 1.44 for ISIS 3He detectors and 1/0.83 for ISIS:LAD circa 1990 scintillator detectors.");
132 auto orderValidator = std::make_shared<Kernel::BoundedValidator<int>>(1, 2);
133 declareProperty(
"Order", 1, orderValidator,
"Placzek correction order (1 or 2), default to 1 (self scattering).");
135 "Sample temperature in Kelvin."
136 "The input property is prioritized over the temperature recorded in the sample log."
137 "The temperature is necessary for computing second order correction.");
138 declareProperty(
"ScaleByPackingFraction",
true,
"Scale the correction value by packing fraction.");
143 "Workspace with the Placzek scattering correction factors.");
153 std::map<std::string, std::string> issues;
159 if (specInfo.
size() == 0) {
160 issues[
"InputWorkspace"] =
"Input workspace does not have detector information";
164 if (
isDefault(
"SampleTemperature") && (order == 2)) {
165 const auto run = inWS->run();
168 if (run.hasProperty(
"SampleTemp")) {
169 sampleTempLog = run.getLogData(
"SampleTemp");
172 if (run.hasProperty(
"sample_temp")) {
173 sampleTempLog = run.getLogData(
"sample_temp");
175 if (!sampleTempLog) {
176 issues[
"SampleTemperature"] =
"Cannot locate sample temperature in the run.";
182 const int64_t numHist = incidentWS->spectrumInfo().size();
186 issues[
"IncidentSpectra"] =
"Need three spectra here for second order calculation.";
189 if (incidentWS->readY(0).empty()) {
190 issues[
"IncidentSpectra"] =
"Flux is empty";
192 if (incidentWS->readY(1).empty()) {
193 issues[
"IncidentSpectra"] =
"First order derivate of the incident spectrum is empty";
195 if (incidentWS->readY(2).empty()) {
196 issues[
"IncidentSpectra"] =
"Second order derivate of the incident spectrum is empty";
202 issues[
"IncidentSpectra"] =
"Need two spectra here for first order calculation.";
205 if (incidentWS->readY(0).empty()) {
206 issues[
"IncidentSpectra"] =
"Flux is empty";
208 if (incidentWS->readY(1).empty()) {
209 issues[
"IncidentSpectra"] =
"First order derivate of the incident spectrum is empty";
217 const int64_t numHistEff = efficiencyWS->spectrumInfo().size();
220 if (numHistEff < 3) {
221 issues[
"EfficiencySpectra"] =
"Need three spectra here for second order calculation.";
224 if (efficiencyWS->readY(0).empty()) {
225 issues[
"EfficiencySpectra"] =
"Detector efficiency is empty";
227 if (efficiencyWS->readY(1).empty()) {
228 issues[
"EfficiencySpectra"] =
"First order derivate of the efficiency spectrum is empty";
230 if (efficiencyWS->readY(2).empty()) {
231 issues[
"EfficiencySpectra"] =
"Second order derivate of the efficiency spectrum is empty";
236 if (numHistEff < 2) {
237 issues[
"EfficiencySpectra"] =
"Need two spectra here for first order calculation.";
240 if (efficiencyWS->readY(0).empty()) {
241 issues[
"EfficiencySpectra"] =
"Detector efficiency is empty";
243 if (efficiencyWS->readY(1).empty()) {
244 issues[
"EfficiencySpectra"] =
"First order derivate of the efficiency spectrum is empty";
262 const bool scaleByPackingFraction =
getProperty(
"ScaleByPackingFraction");
269 DataObjects::create<API::HistoWorkspace>(*inWS, incidentWS->getSpectrum(0).binEdges());
270 outputWS->getAxis(0)->unit() = incidentWS->getAxis(0)->unit();
271 outputWS->setDistribution(
true);
272 outputWS->setYUnit(
"");
273 outputWS->setYUnitLabel(
"Counts");
288 const auto xLambda = incidentWS->getSpectrum(0).points();
291 const std::pair<double, double> summationTerms = calculateSummationTerm(inWS->sample().getMaterial());
309 const std::vector<double> phi2 = (order == 2) ?
getFluxCoefficient2() : std::vector<double>();
314 const int64_t numHist = specInfo.
size();
316 for (int64_t specIndex = 0; specIndex < numHist; specIndex++) {
318 auto &
y = outputWS->mutableY(specIndex);
322 if (!specInfo.
isMonitor(specIndex) && !(specInfo.
l2(specIndex) == 0.0)) {
326 double l1 = specInfo.
l1();
337 const double sinThetaBy2 = sin(
twoTheta / 2.0);
338 const double f = l1 / (l1 +
l2);
340 const double kBT = BoltzmannConstant * sampleTemperature;
342 const double sinHalfAngleSq = sinThetaBy2 * sinThetaBy2;
344 for (
size_t xIndex = 0; xIndex < xLambda.size(); xIndex++) {
346 const double term1 = (f - 1.0) * phi1[xIndex];
347 const double term2 = f * (1.0 - eps1[xIndex]);
348 double inelasticPlaczekCorrection = 2.0 * (term1 + term2 - 3) * sinHalfAngleSq * summationTerms.first;
351 const double k = 2 * M_PI / xLambda[xIndex];
352 const double energy = E_mev_toNeutronWavenumberSq * (k * k);
353 const double kBToverE = kBT /
energy;
355 const double bracket_1 = (8 * f - 9) * (f - 1) * phi1[xIndex]
356 - 3 * f * (2 * f - 3) * eps1[xIndex]
357 + 2 * f * (1 - f) * phi1[xIndex] * eps1[xIndex]
358 + (1 - f) * (1 - f) * phi2[xIndex]
359 + f * f * eps2[xIndex]
360 + 3 * (4 * f - 5) * (f - 1);
361 const double P2_part1 = summationTerms.first * (kBToverE / 2.0 + kBToverE * sinHalfAngleSq * bracket_1);
362 const double bracket_2 = (4 * f - 7) * (f - 1) * phi1[xIndex]
363 + f * (7 - 2 * f) * eps1[xIndex]
364 + 2 * f * (1 - f) * phi1[xIndex] * eps1[xIndex]
365 + (1 - f) * (1 - f) * phi2[xIndex]
366 + f * f * eps2[xIndex]
367 + (2 * f * f - 7 * f + 8);
368 const double P2_part2 = 2 * sinHalfAngleSq * summationTerms.second * (1 + sinHalfAngleSq * bracket_2);
370 inelasticPlaczekCorrection += P2_part1 + P2_part2;
374 scaleByPackingFraction ? (1 + inelasticPlaczekCorrection) * packingFraction : inelasticPlaczekCorrection;
377 for (
size_t xIndex = 0; xIndex < xLambda.size(); xIndex++) {
398 const auto &material = ws->sample().getMaterial();
401 double packingFraction = material.packingFraction();
404 const double crystalDensity =
getProperty(
"CrystalDensity");
405 if (crystalDensity > 0.) {
407 packingFraction = material.numberDensity() / crystalDensity;
410 return packingFraction;
419 double sampleTemperature =
getProperty(
"SampleTemperature");
426 if (
isDefault(
"SampleTemperature") && (order == 2)) {
429 const auto run = inWS->run();
431 if (run.hasProperty(
"SampleTemp")) {
432 sampleTemperature = run.getPropertyAsSingleValue(
"SampleTemp");
433 const std::string sampleTempUnit = run.getProperty(
"SampleTemp")->units();
434 if (sampleTempUnit ==
"C") {
435 sampleTemperature = sampleTemperature + 273.15;
438 }
else if (run.hasProperty(
"sample_temp")) {
439 sampleTemperature = run.getPropertyAsSingleValue(
"sample_temp");
440 const std::string sampleTempUnit = run.getProperty(
"sample_temp")->units();
441 if (sampleTempUnit ==
"C") {
442 sampleTemperature = sampleTemperature + 273.15;
446 throw std::runtime_error(
"Sample temperature is not found in the log.");
449 return sampleTemperature;
459 std::vector<double> phi1;
462 const auto xLambda = incidentWS->getSpectrum(0).points();
463 const auto &incident = incidentWS->readY(0);
464 const auto &incidentPrime = incidentWS->readY(1);
466 for (
size_t i = 0; i < xLambda.size(); i++) {
467 phi1.emplace_back(xLambda[i] * incidentPrime[i] / incident[i]);
479 std::vector<double> phi2;
482 const auto xLambda = incidentWS->getSpectrum(0).points();
483 const auto &incident = incidentWS->readY(0);
484 const auto &incidentPrime2 = incidentWS->readY(2);
486 for (
size_t i = 0; i < xLambda.size(); i++) {
487 phi2.emplace_back(xLambda[i] * xLambda[i] * incidentPrime2[i] / incident[i]);
500 std::vector<double> eps1;
505 const auto xLambda = incidentWS->getSpectrum(0).points();
510 std::vector<double> eps = efficiencyWS->readY(0);
511 std::vector<double> epsPrime = efficiencyWS->readY(1);
512 for (
size_t i = 0; i < xLambda.size(); i++) {
513 double lambda = xLambda[i];
514 double k = 2.0 * M_PI /
lambda;
515 double eps_i = (eps[i] + eps[i + 1]) / 2.0;
516 double epsPrime_i = (epsPrime[i] + epsPrime[i + 1]) / 2.0;
517 eps1.emplace_back(k * epsPrime_i / eps_i);
522 for (
auto x : xLambda) {
524 eps1.emplace_back(
x * exp(
x) / (1.0 - exp(
x)));
537 std::vector<double> eps2;
542 const auto xLambda = incidentWS->getSpectrum(0).points();
548 std::vector<double> eps = efficiencyWS->readY(0);
549 std::vector<double> epsPrime2 = efficiencyWS->readY(2);
550 for (
size_t i = 0; i < xLambda.size(); i++) {
551 double lambda = xLambda[i];
552 double k = 2.0 * M_PI /
lambda;
553 double eps_i = (eps[i] + eps[i + 1]) / 2.0;
554 double epsPrime2_i = (epsPrime2[i] + epsPrime2[i + 1]) / 2.0;
555 eps2.emplace_back(k * k * epsPrime2_i / eps_i);
565 for (
auto x : xLambda) {
567 double eps1 =
x * exp(
x) / (1.0 - exp(
x));
568 eps2.emplace_back((-
x - 2.0) * eps1);
#define DECLARE_ALGORITHM(classname)
const std::vector< double > * lambda
#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.
void declareProperty(std::unique_ptr< Kernel::Property > p, const std::string &doc="") override
Add a property to the list of managed properties.
TypedValue getProperty(const std::string &name) const override
Get the value of a property.
bool isDefault(const std::string &name) const
A validator which checks that a workspace has a valid instrument.
A validator which checks that sample has the required properties.
API::SpectrumInfo is an intermediate step towards a SpectrumInfo that is part of Instrument-2....
bool isMonitor(const size_t index) const
Returns true if the detector(s) associated with the spectrum are monitors.
double l2(const size_t index) const
Returns L2 (distance from sample to spectrum).
size_t size() const
Returns the size of the SpectrumInfo, i.e., the number of spectra.
double l1() const
Returns L1 (distance from source to sample).
void getDetectorValues(const Kernel::Unit &inputUnit, const Kernel::Unit &outputUnit, const Kernel::DeltaEMode::Type emode, const bool signedTheta, int64_t wsIndex, Kernel::UnitParametersMap &pmap) const
Get the detector values relevant to unit conversion for a workspace index.
A property class for workspaces.
A validator which checks that the unit of the workspace referred to by a WorkspaceProperty is the exp...
CalculatePlaczek : Placzek 1st&2nd order correction for inelastic scattering.
std::vector< double > getEfficiencyCoefficient2()
compute the second order detector efficiency coefficient vector
std::map< std::string, std::string > validateInputs() override
validate inputs
const std::string category() const override
Algorithm's category for identification.
void exec() override
Execute the algorithm.
const std::string summary() const override
Algorithm's summary for use in the GUI and help.
const std::vector< std::string > seeAlso() const override
Algorithm's see also for use in the GUI and help.
double getSampleTemperature()
query the sample temperature from input property or sample log
std::vector< double > getFluxCoefficient2()
Compute the flux coefficient for the second order, i.e.
double getPackingFraction(const API::MatrixWorkspace_const_sptr &ws)
compute the packing fraction with given crystal density
std::vector< double > getFluxCoefficient1()
Compute the flux coefficient for the first order, i.e.
int version() const override
Algorithm's version for identification.
std::vector< double > getEfficiencyCoefficient1()
Compute the detector efficiency coefficient based on either given efficiency workspace or a vector de...
void init() override
Initialize the algorithm's properties.
IPropertyManager * setProperty(const std::string &name, const T &value)
Templated method to set the value of a PropertyWithValue.
void information(const std::string &msg)
Logs at information level.
Base class for properties.
void initialize(const double &_l1, const int &_emode, const UnitParametersMap ¶ms)
Initialize the unit to perform conversion using singleToTof() and singleFromTof()
Time of flight in microseconds.
std::shared_ptr< const MatrixWorkspace > MatrixWorkspace_const_sptr
shared pointer to the matrix workspace base class (const version)
std::shared_ptr< MatrixWorkspace > MatrixWorkspace_sptr
shared pointer to the matrix workspace base class
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.
static constexpr double BoltzmannConstant
Boltzmann Constant in meV/K Taken from http://physics.nist.gov/cuu/Constants on 10/07/2012.
static constexpr double E_mev_toNeutronWavenumberSq
Transformation coefficient to transform neutron energy into neutron wavevector: K-neutron[m^-10] = sq...
static constexpr double NeutronMass
Mass of the neutron in kg.
static constexpr double AtomicMassUnit
AMU in kg.
constexpr double EMPTY_DBL() noexcept
Returns what we consider an "empty" double within a property.
Generate a tableworkspace to store the calibration results.
Describes the direction (within an algorithm) of a Property.
@ Input
An input workspace.
@ Output
An output workspace.