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