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