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