Line data Source code
1 : #include "crpropa/module/EMPairProduction.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 : #include "kiss/logger.h"
5 :
6 : #include <fstream>
7 : #include <limits>
8 : #include <stdexcept>
9 :
10 :
11 : namespace crpropa {
12 :
13 : static const double mec2 = mass_electron * c_squared;
14 :
15 9 : EMPairProduction::EMPairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit) {
16 9 : setPhotonField(photonField);
17 9 : setThinning(thinning);
18 9 : setLimit(limit);
19 9 : setHaveElectrons(haveElectrons);
20 9 : }
21 :
22 33 : void EMPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
23 33 : this->photonField = photonField;
24 33 : std::string fname = photonField->getFieldName();
25 33 : setDescription("EMPairProduction: " + fname);
26 66 : initRate(getDataPath("EMPairProduction/rate_" + fname + ".txt"));
27 99 : initCumulativeRate(getDataPath("EMPairProduction/cdf_" + fname + ".txt"));
28 33 : }
29 :
30 14 : void EMPairProduction::setHaveElectrons(bool haveElectrons) {
31 14 : this->haveElectrons = haveElectrons;
32 14 : }
33 :
34 9 : void EMPairProduction::setLimit(double limit) {
35 9 : this->limit = limit;
36 9 : }
37 :
38 13 : void EMPairProduction::setThinning(double thinning) {
39 13 : this->thinning = thinning;
40 13 : }
41 :
42 33 : void EMPairProduction::initRate(std::string filename) {
43 33 : std::ifstream infile(filename.c_str());
44 :
45 33 : if (!infile.good())
46 0 : throw std::runtime_error("EMPairProduction: could not open file " + filename);
47 :
48 : // clear previously loaded interaction rates
49 : tabEnergy.clear();
50 : tabRate.clear();
51 :
52 7510 : while (infile.good()) {
53 7477 : if (infile.peek() != '#') {
54 : double a, b;
55 : infile >> a >> b;
56 7345 : if (infile) {
57 7312 : tabEnergy.push_back(pow(10, a) * eV);
58 7312 : tabRate.push_back(b / Mpc);
59 : }
60 : }
61 7477 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
62 : }
63 33 : infile.close();
64 33 : }
65 :
66 33 : void EMPairProduction::initCumulativeRate(std::string filename) {
67 33 : std::ifstream infile(filename.c_str());
68 :
69 33 : if (!infile.good())
70 0 : throw std::runtime_error("EMPairProduction: could not open file " + filename);
71 :
72 : // clear previously loaded tables
73 : tabE.clear();
74 : tabs.clear();
75 : tabCDF.clear();
76 :
77 : // skip header
78 165 : while (infile.peek() == '#')
79 132 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
80 :
81 : // read s values in first line
82 : double a;
83 : infile >> a; // skip first value
84 3663 : while (infile.good() and (infile.peek() != '\n')) {
85 : infile >> a;
86 3630 : tabs.push_back(pow(10, a) * eV * eV);
87 : }
88 :
89 : // read all following lines: E, cdf values
90 7345 : while (infile.good()) {
91 : infile >> a;
92 7345 : if (!infile)
93 : break; // end of file
94 7312 : tabE.push_back(pow(10, a) * eV);
95 : std::vector<double> cdf;
96 811632 : for (int i = 0; i < tabs.size(); i++) {
97 : infile >> a;
98 804320 : cdf.push_back(a / Mpc);
99 : }
100 7312 : tabCDF.push_back(cdf);
101 7312 : }
102 33 : infile.close();
103 33 : }
104 :
105 : // Hold an data array to interpolate the energy distribution on
106 : class PPSecondariesEnergyDistribution {
107 : private:
108 : std::vector<double> tab_s;
109 : std::vector< std::vector<double> > data;
110 : size_t N;
111 :
112 : public:
113 : // differential cross section for pair production for x = Epositron/Egamma, compare Lee 96 arXiv:9604098
114 : double dSigmadE_PPx(double x, double beta) {
115 1000000 : double A = (x / (1. - x) + (1. - x) / x );
116 1000000 : double B = (1. / x + 1. / (1. - x) );
117 1000 : double y = (1 - beta * beta);
118 1000000 : return A + y * B - y * y / 4 * B * B;
119 : }
120 :
121 1 : PPSecondariesEnergyDistribution() {
122 1 : N = 1000;
123 : size_t Ns = 1000;
124 : double s_min = 4 * mec2 * mec2;
125 : double s_max = 1e23 * eV * eV;
126 : double dls = log(s_max / s_min) / Ns;
127 1 : data = std::vector< std::vector<double> >(Ns, std::vector<double>(N));
128 1 : tab_s = std::vector<double>(Ns + 1);
129 :
130 1002 : for (size_t i = 0; i < Ns + 1; ++i)
131 1001 : tab_s[i] = s_min * exp(i*dls); // tabulate s bin borders
132 :
133 1001 : for (size_t i = 0; i < Ns; i++) {
134 1000 : double s = s_min * exp(i*dls + 0.5*dls);
135 1000 : double beta = sqrt(1 - s_min/s);
136 1000 : double x0 = (1 - beta) / 2;
137 1000 : double dx = log((1 + beta) / (1 - beta)) / N;
138 :
139 : // cumulative midpoint integration
140 1000 : std::vector<double> data_i(1000);
141 1000 : data_i[0] = dSigmadE_PPx(x0, beta) * expm1(dx);
142 1000000 : for (size_t j = 1; j < N; j++) {
143 999000 : double x = x0 * exp(j*dx + 0.5*dx);
144 999000 : double binWidth = exp((j+1)*dx)-exp(j*dx);
145 999000 : data_i[j] = dSigmadE_PPx(x, beta) * binWidth + data_i[j-1];
146 : }
147 1000 : data[i] = data_i;
148 1000 : }
149 1 : }
150 :
151 : // sample positron energy from cdf(E, s_kin)
152 561 : double sample(double E0, double s) {
153 : // get distribution for given s
154 561 : size_t idx = std::lower_bound(tab_s.begin(), tab_s.end(), s) - tab_s.begin();
155 561 : if (idx > data.size())
156 : return NAN;
157 :
158 561 : std::vector<double> s0 = data[idx];
159 :
160 : // draw random bin
161 561 : Random &random = Random::instance();
162 561 : size_t j = random.randBin(s0) + 1;
163 :
164 : double s_min = 4. * mec2 * mec2;
165 561 : double beta = sqrtl(1. - s_min / s);
166 561 : double x0 = (1. - beta) / 2.;
167 561 : double dx = log((1 + beta) / (1 - beta)) / N;
168 561 : double binWidth = x0 * (exp(j*dx) - exp((j-1)*dx));
169 561 : if (random.rand() < 0.5)
170 285 : return E0 * (x0 * exp((j-1) * dx) + binWidth);
171 : else
172 276 : return E0 * (1 - (x0 * exp((j-1) * dx) + binWidth));
173 561 : }
174 : };
175 :
176 561 : void EMPairProduction::performInteraction(Candidate *candidate) const {
177 :
178 : // photon is lost after interacting
179 561 : candidate->setActive(false);
180 :
181 : // check if secondary electron pair needs to be produced
182 561 : if (not haveElectrons)
183 0 : return;
184 :
185 : // scale particle energy instead of background photon energy
186 561 : double z = candidate->getRedshift();
187 561 : double E = candidate->current.getEnergy() * (1 + z);
188 :
189 : // check if in tabulated energy range
190 561 : if (E < tabE.front() or (E > tabE.back())) {
191 0 : KISS_LOG_WARNING
192 0 : << "EMPairProduction: Energy "
193 0 : << E / eV << " eV is not in tabulated range";
194 0 : return;
195 : }
196 :
197 : // sample the value of s
198 561 : Random &random = Random::instance();
199 561 : size_t i = closestIndex(E, tabE); // find closest tabulation point
200 561 : size_t j = random.randBin(tabCDF[i]);
201 561 : if (j <= 0) {
202 84 : KISS_LOG_WARNING
203 42 : << "EMPaiProduction: Sampled s value is the lowest tabulated value, which is not physical."
204 42 : << " The index j will be set to 1 to avoid division by zero.";
205 : j = 1; // ensure j is at least 1 to avoid division by
206 : }
207 1122 : double lo = std::max(4 * mec2 * mec2, tabs[j-1]); // first s-tabulation point below min(s_kin) = (2 me c^2)^2; ensure physical value
208 561 : double hi = tabs[j];
209 561 : double s = lo + random.rand() * (hi - lo);
210 :
211 : // sample electron / positron energy
212 561 : static PPSecondariesEnergyDistribution interpolation;
213 561 : double Ee = interpolation.sample(E, s);
214 561 : double Ep = E - Ee;
215 561 : double f = Ep / E;
216 :
217 : // for some backgrounds Ee=nan due to precision limitations.
218 561 : if (not std::isfinite(Ee) || not std::isfinite(Ep)) {
219 0 : KISS_LOG_WARNING
220 0 : << "EMPairProduction: Sampled energies are not finite for primary energy "
221 0 : << E / eV << " eV and s = " << s / (eV * eV) << " eV^2 (maximum tabulated s = "
222 0 : << tabs.back() / (eV * eV) << " eV^2).";
223 0 : return;
224 : }
225 :
226 : // sample random position along current step
227 561 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
228 : // apply sampling
229 561 : if (random.rand() < pow(f, thinning)) {
230 561 : double w = 1. / pow(f, thinning);
231 561 : candidate->addSecondary(11, Ep / (1 + z), pos, w, interactionTag);
232 : }
233 561 : if (random.rand() < pow(1 - f, thinning)){
234 561 : double w = 1. / pow(1 - f, thinning);
235 561 : candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag);
236 : }
237 : }
238 :
239 2549 : void EMPairProduction::process(Candidate *candidate) const {
240 : // check if photon
241 2549 : if (candidate->current.getId() != 22)
242 : return;
243 :
244 : // scale particle energy instead of background photon energy
245 841 : double z = candidate->getRedshift();
246 841 : double E = candidate->current.getEnergy() * (1 + z);
247 :
248 : // check if in tabulated energy range
249 841 : if ((E < tabEnergy.front()) or (E > tabEnergy.back()))
250 : return;
251 :
252 : // interaction rate
253 651 : double rate = interpolate(E, tabEnergy, tabRate);
254 651 : rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
255 :
256 : // run this loop at least once to limit the step size
257 651 : double step = candidate->getCurrentStep();
258 651 : Random &random = Random::instance();
259 : do {
260 651 : double randDistance = -log(random.rand()) / rate;
261 : // check for interaction; if it doesn't ocurr, limit next step
262 651 : if (step < randDistance) {
263 91 : candidate->limitNextStep(limit / rate);
264 : } else {
265 560 : performInteraction(candidate);
266 560 : return;
267 : }
268 91 : step -= randDistance;
269 91 : } while (step > 0.);
270 :
271 : }
272 :
273 1 : void EMPairProduction::setInteractionTag(std::string tag) {
274 1 : interactionTag = tag;
275 1 : }
276 :
277 2 : std::string EMPairProduction::getInteractionTag() const {
278 2 : return interactionTag;
279 : }
280 :
281 : } // namespace crpropa
|
