Line data Source code
1 : #include "crpropa/module/PhotoDisintegration.h"
2 : #include "crpropa/Common.h"
3 : #include "crpropa/Units.h"
4 : #include "crpropa/ParticleID.h"
5 : #include "crpropa/ParticleMass.h"
6 : #include "crpropa/Random.h"
7 : #include "kiss/logger.h"
8 :
9 : #include <cmath>
10 : #include <limits>
11 : #include <sstream>
12 : #include <fstream>
13 : #include <stdexcept>
14 :
15 : namespace crpropa {
16 :
17 : const double PhotoDisintegration::lgmin = 6; // minimum log10(Lorentz-factor)
18 : const double PhotoDisintegration::lgmax = 14; // maximum log10(Lorentz-factor)
19 : const size_t PhotoDisintegration::nlg = 201; // number of Lorentz-factor steps
20 :
21 21 : PhotoDisintegration::PhotoDisintegration(ref_ptr<PhotonField> f, bool havePhotons, double limit, bool superheavy) {
22 42 : setPhotonField(f, superheavy);
23 21 : this->havePhotons = havePhotons;
24 21 : this->limit = limit;
25 21 : }
26 :
27 34 : void PhotoDisintegration::setPhotonField(ref_ptr<PhotonField> photonField, bool superheavy) {
28 34 : this->photonField = photonField;
29 34 : std::string fname = photonField->getFieldName();
30 34 : setDescription("PhotoDisintegration: " + fname);
31 34 : if (superheavy) {
32 10 : initRate(getDataPath("Photodisintegration/rate_" + fname + "_superheavy.txt"));
33 15 : initBranching(getDataPath("Photodisintegration/branching_" + fname + "_superheavy.txt"));
34 : } else {
35 58 : initRate(getDataPath("Photodisintegration/rate_" + fname + ".txt"));
36 87 : initBranching(getDataPath("Photodisintegration/branching_" + fname + ".txt"));
37 : }
38 102 : initPhotonEmission(getDataPath("Photodisintegration/photon_emission_" + fname.substr(0,3) + ".txt"));
39 34 : }
40 :
41 1 : void PhotoDisintegration::setHavePhotons(bool havePhotons) {
42 1 : this->havePhotons = havePhotons;
43 1 : }
44 :
45 0 : void PhotoDisintegration::setLimit(double limit) {
46 0 : this->limit = limit;
47 0 : }
48 :
49 34 : void PhotoDisintegration::initRate(std::string filename) {
50 34 : std::ifstream infile(filename.c_str());
51 34 : if (not infile.good())
52 0 : throw std::runtime_error("PhotoDisintegration: could not open file " + filename);
53 :
54 : pdRate.clear();
55 34 : pdRate.resize((NUCLEAR_ZMAX + 1) * NUCLEAR_NSTRIDE);
56 :
57 : std::string line;
58 11451 : while (std::getline(infile, line)) {
59 11417 : if (line[0] == '#')
60 102 : continue;
61 11315 : std::stringstream lineStream(line);
62 :
63 : int Z, N;
64 11315 : lineStream >> Z >> N;
65 :
66 : double r;
67 2285630 : for (size_t i = 0; i < nlg; i++) {
68 : lineStream >> r;
69 2274315 : pdRate[Z * NUCLEAR_NSTRIDE + N].push_back(r / Mpc);
70 : }
71 11315 : }
72 34 : infile.close();
73 34 : }
74 :
75 34 : void PhotoDisintegration::initBranching(std::string filename) {
76 34 : std::ifstream infile(filename.c_str());
77 34 : if (not infile.good())
78 0 : throw std::runtime_error("PhotoDisintegration: could not open file " + filename);
79 :
80 : pdBranch.clear();
81 34 : pdBranch.resize((NUCLEAR_ZMAX + 1) * NUCLEAR_NSTRIDE);
82 :
83 : std::string line;
84 143879 : while (std::getline(infile, line)) {
85 143845 : if (line[0] == '#')
86 102 : continue;
87 :
88 143743 : std::stringstream lineStream(line);
89 :
90 : int Z, N;
91 143743 : lineStream >> Z >> N;
92 :
93 : Branch branch;
94 143743 : lineStream >> branch.channel;
95 :
96 : double r;
97 29036086 : for (size_t i = 0; i < nlg; i++) {
98 : lineStream >> r;
99 28892343 : branch.branchingRatio.push_back(r);
100 : }
101 :
102 143743 : pdBranch[Z * NUCLEAR_NSTRIDE + N].push_back(branch);
103 143743 : }
104 :
105 34 : infile.close();
106 34 : }
107 :
108 34 : void PhotoDisintegration::initPhotonEmission(std::string filename) {
109 34 : std::ifstream infile(filename.c_str());
110 34 : if (not infile.good())
111 0 : throw std::runtime_error("PhotoDisintegration: could not open file " + filename);
112 :
113 : // clear previously loaded emission probabilities
114 : pdPhoton.clear();
115 :
116 : std::string line;
117 323238 : while (std::getline(infile, line)) {
118 323204 : if (line[0] == '#')
119 102 : continue;
120 :
121 323102 : std::stringstream lineStream(line);
122 :
123 : int Z, N, Zd, Nd;
124 323102 : lineStream >> Z;
125 323102 : lineStream >> N;
126 323102 : lineStream >> Zd;
127 323102 : lineStream >> Nd;
128 :
129 : PhotonEmission em;
130 : lineStream >> em.energy;
131 323102 : em.energy *= eV;
132 :
133 : double r;
134 65266604 : for (size_t i = 0; i < nlg; i++) {
135 : lineStream >> r;
136 64943502 : em.emissionProbability.push_back(r);
137 : }
138 :
139 323102 : int key = Z * 1000000 + N * 10000 + Zd * 100 + Nd;
140 323102 : if (pdPhoton.find(key) == pdPhoton.end()) {
141 : std::vector<PhotonEmission> emissions;
142 63512 : pdPhoton[key] = emissions;
143 63512 : }
144 323102 : pdPhoton[key].push_back(em);
145 323102 : }
146 :
147 34 : infile.close();
148 34 : }
149 :
150 4558 : void PhotoDisintegration::process(Candidate *candidate) const {
151 : // execute the loop at least once for limiting the next step
152 4558 : double step = candidate->getCurrentStep();
153 : do {
154 : // check if nucleus
155 4804 : int id = candidate->current.getId();
156 4804 : if (not isNucleus(id))
157 : return;
158 :
159 4803 : int A = massNumber(id);
160 4803 : int Z = chargeNumber(id);
161 4803 : int N = A - Z;
162 4803 : size_t idx = Z * NUCLEAR_NSTRIDE + N;
163 :
164 : // check if disintegration data available
165 4803 : if ((Z > NUCLEAR_ZMAX) or (N > NUCLEAR_NMAX))
166 : return;
167 4801 : if (pdRate[idx].size() == 0)
168 : return;
169 :
170 : // check if in tabulated energy range
171 602 : double z = candidate->getRedshift();
172 602 : double lg = log10(candidate->current.getLorentzFactor() * (1 + z));
173 602 : if ((lg <= lgmin) or (lg >= lgmax))
174 : return;
175 :
176 602 : double rate = interpolateEquidistant(lg, lgmin, lgmax, pdRate[idx]);
177 602 : rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z); // cosmological scaling, rate per comoving distance
178 :
179 : // check if interaction occurs in this step
180 : // otherwise limit next step to a fraction of the mean free path
181 602 : Random &random = Random::instance();
182 602 : double randDist = -log(random.rand()) / rate;
183 602 : if (step < randDist) {
184 356 : candidate->limitNextStep(limit / rate);
185 356 : return;
186 : }
187 :
188 : // select channel and interact
189 : const std::vector<Branch> &branches = pdBranch[idx];
190 246 : double cmp = random.rand();
191 246 : int l = round((lg - lgmin) / (lgmax - lgmin) * (nlg - 1)); // index of closest tabulation point
192 : size_t i = 0;
193 2015 : while ((i < branches.size()) and (cmp > 0)) {
194 1769 : cmp -= branches[i].branchingRatio[l];
195 1769 : i++;
196 : }
197 246 : performInteraction(candidate, branches[i-1].channel);
198 :
199 : // repeat with remaining step
200 246 : step -= randDist;
201 246 : } while (step > 0);
202 : }
203 :
204 248 : void PhotoDisintegration::performInteraction(Candidate *candidate, int channel) const {
205 248 : KISS_LOG_DEBUG << "Photodisintegration::performInteraction. Channel " << channel << " on candidate " << candidate->getDescription();
206 : // parse disintegration channel
207 : int nNeutron = digit(channel, 100000);
208 : int nProton = digit(channel, 10000);
209 : int nH2 = digit(channel, 1000);
210 : int nH3 = digit(channel, 100);
211 : int nHe3 = digit(channel, 10);
212 : int nHe4 = digit(channel, 1);
213 :
214 248 : int dA = -nNeutron - nProton - 2 * nH2 - 3 * nH3 - 3 * nHe3 - 4 * nHe4;
215 248 : int dZ = -nProton - nH2 - nH3 - 2 * nHe3 - 2 * nHe4;
216 :
217 248 : int id = candidate->current.getId();
218 248 : int A = massNumber(id);
219 248 : int Z = chargeNumber(id);
220 248 : double EpA = candidate->current.getEnergy() / A;
221 :
222 : // create secondaries
223 248 : Random &random = Random::instance();
224 248 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
225 : try
226 : {
227 452 : for (size_t i = 0; i < nNeutron; i++)
228 204 : candidate->addSecondary(nucleusId(1, 0), EpA, pos, 1., interactionTag);
229 345 : for (size_t i = 0; i < nProton; i++)
230 97 : candidate->addSecondary(nucleusId(1, 1), EpA, pos, 1., interactionTag);
231 248 : for (size_t i = 0; i < nH2; i++)
232 0 : candidate->addSecondary(nucleusId(2, 1), EpA * 2, pos, 1., interactionTag);
233 250 : for (size_t i = 0; i < nH3; i++)
234 2 : candidate->addSecondary(nucleusId(3, 1), EpA * 3, pos, 1., interactionTag);
235 248 : for (size_t i = 0; i < nHe3; i++)
236 0 : candidate->addSecondary(nucleusId(3, 2), EpA * 3, pos, 1., interactionTag);
237 289 : for (size_t i = 0; i < nHe4; i++)
238 41 : candidate->addSecondary(nucleusId(4, 2), EpA * 4, pos, 1., interactionTag);
239 :
240 :
241 : // update particle
242 : candidate->created = candidate->current;
243 248 : candidate->current.setId(nucleusId(A + dA, Z + dZ));
244 248 : candidate->current.setEnergy(EpA * (A + dA));
245 : }
246 0 : catch (std::runtime_error &e)
247 : {
248 0 : KISS_LOG_ERROR << "Something went wrong in the PhotoDisentigration\n" << "Please report this error on https://github.com/CRPropa/CRPropa3/issues including your simulation setup and the following random seed:\n" << Random::instance().getSeed_base64();
249 0 : throw;
250 0 : }
251 :
252 248 : if (not havePhotons)
253 : return;
254 :
255 : // create photons
256 29 : double z = candidate->getRedshift();
257 29 : double lg = log10(candidate->current.getLorentzFactor() * (1 + z));
258 29 : double lf = candidate->current.getLorentzFactor();
259 :
260 29 : int l = round((lg - lgmin) / (lgmax - lgmin) * (nlg - 1)); // index of closest tabulation point
261 29 : int key = Z*1e6 + (A-Z)*1e4 + (Z+dZ)*1e2 + (A+dA) - (Z+dZ);
262 :
263 185 : for (int i = 0; i < pdPhoton[key].size(); i++) {
264 : // check for random emission
265 156 : if (random.rand() > pdPhoton[key][i].emissionProbability[l])
266 131 : continue;
267 :
268 : // boost to lab frame
269 25 : double cosTheta = 2 * random.rand() - 1;
270 25 : double E = pdPhoton[key][i].energy * lf * (1 - cosTheta);
271 25 : candidate->addSecondary(22, E, pos, 1., interactionTag);
272 : }
273 : }
274 :
275 3528 : double PhotoDisintegration::lossLength(int id, double gamma, double z) {
276 : // check if nucleus
277 3528 : if (not (isNucleus(id)))
278 : return std::numeric_limits<double>::max();
279 :
280 3528 : int A = massNumber(id);
281 3528 : int Z = chargeNumber(id);
282 3528 : int N = A - Z;
283 3528 : size_t idx = Z * NUCLEAR_NSTRIDE + N;
284 :
285 : // check if disintegration data available
286 3528 : if ((Z > NUCLEAR_ZMAX) or (N > NUCLEAR_NMAX))
287 : return std::numeric_limits<double>::max();
288 : const std::vector<double> &rate = pdRate[idx];
289 3528 : if (rate.size() == 0)
290 : return std::numeric_limits<double>::max();
291 :
292 : // check if in tabulated energy range
293 3528 : double lg = log10(gamma * (1 + z));
294 3528 : if ((lg <= lgmin) or (lg >= lgmax))
295 : return std::numeric_limits<double>::max();
296 :
297 : // total interaction rate
298 3500 : double lossRate = interpolateEquidistant(lg, lgmin, lgmax, rate);
299 :
300 : // comological scaling, rate per physical distance
301 3500 : lossRate *= pow_integer<3>(1 + z) * photonField->getRedshiftScaling(z);
302 :
303 : // average number of nucleons lost for all disintegration channels
304 : double avg_dA = 0;
305 : const std::vector<Branch> &branches = pdBranch[idx];
306 47250 : for (size_t i = 0; i < branches.size(); i++) {
307 43750 : int channel = branches[i].channel;
308 : int dA = 0;
309 : dA += 1 * digit(channel, 100000);
310 43750 : dA += 1 * digit(channel, 10000);
311 43750 : dA += 2 * digit(channel, 1000);
312 43750 : dA += 3 * digit(channel, 100);
313 43750 : dA += 3 * digit(channel, 10);
314 43750 : dA += 4 * digit(channel, 1);
315 :
316 43750 : double br = interpolateEquidistant(lg, lgmin, lgmax, branches[i].branchingRatio);
317 43750 : avg_dA += br * dA;
318 : }
319 :
320 3500 : lossRate *= avg_dA / A;
321 3500 : return 1 / lossRate;
322 : }
323 :
324 1 : void PhotoDisintegration::setInteractionTag(std::string tag) {
325 1 : interactionTag = tag;
326 1 : }
327 :
328 2 : std::string PhotoDisintegration::getInteractionTag() const {
329 2 : return interactionTag;
330 : }
331 :
332 : } // namespace crpropa
|
