1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
|
// ************************************************************************************************
//
// BornAgain: simulate and fit reflection and scattering
//
//! @file Sim/Simulation/ScatteringSimulation.cpp
//! @brief Implements interface ISimulation.
//!
//! @homepage http://www.bornagainproject.org
//! @license GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
// ************************************************************************************************
#include "Sim/Simulation/ScatteringSimulation.h"
#include "Base/Axis/Frame.h"
#include "Base/Axis/Pixel.h"
#include "Base/Axis/Scale.h"
#include "Base/Progress/ProgressHandler.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Beam/IFootprint.h"
#include "Device/Data/Datafield.h"
#include "Device/Detector/IDetector.h"
#include "Param/Distrib/DistributionHandler.h"
#include "Resample/Element/DiffuseElement.h"
#include "Sample/Multilayer/Sample.h"
#include "Sim/Background/IBackground.h"
#include "Sim/Computation/DWBAComputation.h"
ScatteringSimulation::ScatteringSimulation(const Beam& beam, const Sample& sample,
const IDetector& detector)
: ISimulation(sample)
, m_beam(beam.clone())
, m_detector(detector.clone())
{
}
ScatteringSimulation::~ScatteringSimulation() = default;
void ScatteringSimulation::addParameterDistribution(
ParameterDistribution::WhichParameter whichParameter, const IDistribution1D& distribution)
{
ParameterDistribution par_distr(whichParameter, distribution);
distributionHandler().addDistribution(par_distr);
}
//... Overridden executors:
//! init callbacks for setting the parameter values
void ScatteringSimulation::initDistributionHandler()
{
for (const auto& distribution : distributionHandler().paramDistributions()) {
switch (distribution.whichParameter()) {
case ParameterDistribution::BeamAzimuthalAngle:
distributionHandler().defineCallbackForDistribution(
&distribution, [&](double d) { m_beam->setAzimuthalAngle(d); });
break;
case ParameterDistribution::BeamGrazingAngle:
distributionHandler().defineCallbackForDistribution(
&distribution, [&](double d) { m_beam->setGrazingAngle(d); });
break;
case ParameterDistribution::BeamWavelength:
distributionHandler().defineCallbackForDistribution(
&distribution, [&](double d) { m_beam->setWavelength(d); });
break;
default:
ASSERT_NEVER;
}
}
}
void ScatteringSimulation::prepareSimulation()
{
sample()->checkMaterials(beam().wavelength());
m_active_indices = m_detector->activeIndices();
m_pixels.reserve(m_active_indices.size());
for (size_t i : m_active_indices)
m_pixels.push_back(m_detector->createPixel(m_detector->roiToFullIndex(i)));
}
void ScatteringSimulation::runComputation(const ReSample& re_sample, size_t i, double weight)
{
if (m_cache.empty())
m_cache.resize(nElements(), 0.0);
const bool isSpecular = m_active_indices[i] == m_detector->indexOfSpecular(beam());
DiffuseElement ele(beam().wavelength(), beam().alpha_i(), beam().phi_i(), m_pixels[i],
beam().polMatrix(), m_detector->analyzer().matrix(), isSpecular);
double intensity = Compute::scattered_and_reflected(re_sample, options(), ele);
if (const auto* footprint = beam().footprint())
intensity *= footprint->calculate(beam().alpha_i());
double sin_alpha_i = std::abs(std::sin(beam().alpha_i()));
if (sin_alpha_i == 0.0) {
intensity = 0;
} else {
const double solid_angle = m_pixels[i]->solidAngle();
intensity *= m_beam->intensity() * solid_angle / sin_alpha_i;
}
m_cache[i] += intensity * weight;
progress().incrementDone(1);
}
//... Overridden getters:
bool ScatteringSimulation::force_polarized() const
{
return m_detector->analyzer().BlochVector() != R3{};
}
size_t ScatteringSimulation::nElements() const
{
return m_active_indices.size();
}
Datafield ScatteringSimulation::packResult() const
{
Datafield result(m_detector->createDetectorMap());
for (size_t i = 0; i < m_active_indices.size(); ++i)
result.setAt(m_active_indices[i], m_cache[i]);
m_detector->applyDetectorResolution(&result);
if (background())
for (size_t i : m_active_indices)
result[i] = background()->addBackground(result[i]);
return result;
}
|
