Line data Source code
1 : #include "crpropa/module/EMInverseComptonScattering.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 : #include "crpropa/Common.h"
5 :
6 : #include <vector>
7 : #include <cmath>
8 :
9 : namespace crpropa {
10 :
11 : static const double mec2 = mass_electron * c_squared;
12 :
13 5 : EMInverseComptonScattering::EMInverseComptonScattering(ref_ptr<PhotonField> photonField, bool havePhotons, double thinning, double limit, ref_ptr<Surface> surface) {
14 10 : setSurface(surface);
15 5 : setPhotonField(photonField);
16 5 : setHavePhotons(havePhotons);
17 5 : setLimit(limit);
18 5 : setThinning(thinning);
19 5 : }
20 :
21 17 : void EMInverseComptonScattering::setPhotonField(ref_ptr<PhotonField> photonField) {
22 :
23 17 : this->photonField = photonField;
24 17 : std::string fname = photonField->getFieldName();
25 34 : setDescription("EMInverseComptonScattering: " + fname);
26 :
27 17 : if (!this->photonField->hasPositionDependence()) {
28 : this->interactionRates = new InteractionRatesHomogeneous(
29 51 : getDataPath("EMInverseComptonScattering/rate_" + fname + ".txt"),
30 51 : getDataPath("EMInverseComptonScattering/cdf_" + fname + ".txt")
31 34 : );
32 : } else {
33 : this->interactionRates = new InteractionRatesPositionDependent(
34 0 : getDataPath("EMInverseComptonScattering/" + fname + "/Rate/"),
35 0 : getDataPath("EMInverseComptonScattering/" + fname + "/CumulativeRate/"),
36 : this->surface
37 0 : );
38 : }
39 17 : }
40 :
41 6 : void EMInverseComptonScattering::setHavePhotons(bool havePhotons) {
42 6 : this->havePhotons = havePhotons;
43 6 : }
44 :
45 5 : void EMInverseComptonScattering::setLimit(double limit) {
46 5 : this->limit = limit;
47 5 : }
48 :
49 5 : void EMInverseComptonScattering::setThinning(double thinning) {
50 5 : this->thinning = thinning;
51 5 : }
52 :
53 5 : void EMInverseComptonScattering::setSurface(ref_ptr<Surface> surface) {
54 5 : this->surface = surface;
55 5 : }
56 :
57 0 : ref_ptr<Surface> EMInverseComptonScattering::getSurface() const {
58 0 : return this->surface;
59 : }
60 :
61 0 : void EMInverseComptonScattering::setInteractionRates(ref_ptr<InteractionRates> intRates) {
62 0 : this->interactionRates = intRates;
63 0 : }
64 :
65 0 : ref_ptr<InteractionRates> EMInverseComptonScattering::getInteractionRates() const {
66 0 : return this->interactionRates;
67 : }
68 :
69 0 : void EMInverseComptonScattering::initRate(std::string path) {
70 0 : this->interactionRates->initRate(path);
71 0 : }
72 :
73 0 : void EMInverseComptonScattering::initCumulativeRate(std::string path) {
74 0 : this->interactionRates->initCumulativeRate(path);
75 0 : }
76 :
77 : // Class to calculate the energy distribution of the ICS photon and to sample from it
78 : class ICSSecondariesEnergyDistribution {
79 : private:
80 : std::vector< std::vector<double> > data;
81 : std::vector<double> s_values;
82 : size_t Ns;
83 : size_t Nrer;
84 : double s_min;
85 : double s_max;
86 : double dls;
87 :
88 : public:
89 : // differential cross-section, see Lee '96 (arXiv:9604098), eq. 23 for x = Ee'/Ee
90 : double dSigmadE(double x, double beta) {
91 1000000 : double q = ((1 - beta) / beta) * (1 - 1./x);
92 1000000 : return ((1 + beta) / beta) * (x + 1./x + 2 * q + q * q);
93 : }
94 :
95 : // create the cumulative energy distribution of the up-scattered photon
96 1 : ICSSecondariesEnergyDistribution() {
97 1 : Ns = 1000;
98 1 : Nrer = 1000;
99 1 : s_min = mec2 * mec2;
100 1 : s_max = 2e23 * eV * eV;
101 1 : dls = (log(s_max) - log(s_min)) / Ns;
102 1 : data = std::vector< std::vector<double> >(1000, std::vector<double>(1000));
103 1 : std::vector<double> data_i(1000);
104 :
105 : // tabulate s bin borders
106 1 : s_values = std::vector<double>(1001);
107 1002 : for (size_t i = 0; i < Ns + 1; ++i)
108 1001 : s_values[i] = s_min * exp(i*dls);
109 :
110 :
111 : // for each s tabulate cumulative differential cross section
112 1001 : for (size_t i = 0; i < Ns; i++) {
113 1000 : double s = s_min * exp((i+0.5) * dls);
114 1000 : double beta = (s - s_min) / (s + s_min);
115 1000 : double x0 = (1 - beta) / (1 + beta);
116 1000 : double dlx = -log(x0) / Nrer;
117 :
118 : // cumulative midpoint integration
119 1000 : data_i[0] = dSigmadE(x0, beta) * expm1(dlx);
120 1000000 : for (size_t j = 1; j < Nrer; j++) {
121 999000 : double x = x0 * exp((j+0.5) * dlx);
122 999000 : double dx = exp((j+1) * dlx) - exp(j * dlx);
123 999000 : data_i[j] = dSigmadE(x, beta) * dx;
124 999000 : data_i[j] += data_i[j-1];
125 : }
126 1000 : data[i] = data_i;
127 : }
128 1 : }
129 :
130 : // draw random energy for the up-scattered photon Ep(Ee, s)
131 1 : double sample(double Ee, double s) {
132 1 : size_t idx = std::lower_bound(s_values.begin(), s_values.end(), s) - s_values.begin();
133 1 : std::vector<double> s0 = data[idx];
134 1 : Random &random = Random::instance();
135 1 : size_t j = random.randBin(s0) + 1; // draw random bin (upper bin boundary returned)
136 1 : double beta = (s - s_min) / (s + s_min);
137 1 : double x0 = (1 - beta) / (1 + beta);
138 1 : double dlx = -log(x0) / Nrer;
139 1 : double binWidth = x0 * (exp(j * dlx) - exp((j-1) * dlx));
140 1 : double Ep = (x0 * exp((j-1) * dlx) + binWidth) * Ee;
141 1 : return std::min(Ee, Ep); // prevent Ep > Ee from numerical inaccuracies
142 1 : }
143 : };
144 :
145 1 : double EMInverseComptonScattering::getRate(double E, const Vector3d &position, double z) const {
146 1 : return this->interactionRates->getProcessRate(E, position) * pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
147 : }
148 :
149 1 : void EMInverseComptonScattering::performInteraction(Candidate *candidate) const {
150 :
151 : // scale the particle energy instead of background photons
152 1 : double z = candidate->getRedshift();
153 1 : double E = candidate->current.getEnergy() * (1 + z);
154 :
155 1 : Vector3d position = candidate->current.getPosition();
156 :
157 : std::vector<double> tabE;
158 : std::vector<double> tabs;
159 : std::vector<std::vector<double>> tabCDF;
160 :
161 1 : this->interactionRates->loadPerformInteractionTabs(position, tabE, tabs, tabCDF);
162 :
163 1 : if (E < tabE.front() or E > tabE.back())
164 : return;
165 :
166 : // sample the value of s
167 1 : Random &random = Random::instance();
168 1 : size_t i = closestIndex(E, tabE);
169 1 : size_t j = random.randBin(tabCDF[i]);
170 1 : double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1);
171 1 : double s = s_kin + mec2 * mec2;
172 :
173 : // sample electron energy after scattering
174 1 : static ICSSecondariesEnergyDistribution distribution;
175 1 : double Enew = distribution.sample(E, s);
176 :
177 : // add up-scattered photon
178 1 : double Esecondary = E - Enew;
179 1 : double f = Enew / E;
180 1 : if (havePhotons) {
181 1 : if (random.rand() < pow(1 - f, thinning)) {
182 1 : double w = 1. / pow(1 - f, thinning);
183 1 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
184 1 : candidate->addSecondary(22, Esecondary / (1 + z), pos, w, interactionTag);
185 : }
186 : }
187 :
188 : // update the primary particle energy; do this after adding the secondary to correctly set the secondary's parent
189 1 : candidate->current.setEnergy(Enew / (1 + z));
190 1 : }
191 :
192 1307 : void EMInverseComptonScattering::process(Candidate *candidate) const {
193 : // check if electron / positron
194 1307 : int id = candidate->current.getId();
195 1307 : if (abs(id) != 11)
196 1307 : return;
197 :
198 : // scale the particle energy instead of background photons
199 1 : double z = candidate->getRedshift();
200 1 : double E = candidate->current.getEnergy() * (1 + z);
201 1 : Vector3d position = candidate->current.getPosition();
202 :
203 : // interaction rate
204 1 : double rate = getRate(E, position, z);
205 1 : if (rate < 0)
206 : return;
207 :
208 1 : rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
209 :
210 : // run this loop at least once to limit the step size
211 1 : double step = candidate->getCurrentStep();
212 1 : Random &random = Random::instance();
213 : do {
214 1 : double randDistance = -log(random.rand()) / rate;
215 :
216 : // check for interaction; if it doesn't ocurr, limit next step
217 1 : if (step < randDistance) {
218 1 : candidate->limitNextStep(limit / rate);
219 1 : return;
220 : }
221 0 : performInteraction(candidate);
222 :
223 : // repeat with remaining step
224 0 : step -= randDistance;
225 0 : } while (step > 0);
226 : }
227 :
228 1 : void EMInverseComptonScattering::setInteractionTag(std::string tag) {
229 1 : interactionTag = tag;
230 1 : }
231 :
232 2 : std::string EMInverseComptonScattering::getInteractionTag() const {
233 2 : return interactionTag;
234 : }
235 :
236 : } // namespace crpropa
|
