Line data Source code
1 : #include "crpropa/module/NuclearDecay.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/ParticleID.h"
4 : #include "crpropa/ParticleMass.h"
5 : #include "crpropa/Random.h"
6 :
7 : #include <fstream>
8 : #include <limits>
9 : #include <cmath>
10 : #include <stdexcept>
11 :
12 : #include "kiss/logger.h"
13 :
14 : namespace crpropa {
15 :
16 9 : NuclearDecay::NuclearDecay(bool electrons, bool photons, bool neutrinos, double l) {
17 9 : haveElectrons = electrons;
18 9 : havePhotons = photons;
19 9 : haveNeutrinos = neutrinos;
20 9 : limit = l;
21 9 : setDescription("NuclearDecay");
22 :
23 : // load decay table
24 18 : std::string filename = getDataPath("nuclear_decay.txt");
25 9 : std::ifstream infile(filename.c_str());
26 9 : if (!infile.good())
27 0 : throw std::runtime_error(
28 0 : "crpropa::NuclearDecay: could not open file " + filename);
29 :
30 9 : decayTable.resize(27 * 31);
31 : std::string line;
32 8550 : while (std::getline(infile,line)) {
33 8541 : std::stringstream stream(line);
34 8541 : if (stream.peek() == '#')
35 : continue;
36 : DecayMode decay;
37 : int Z, N;
38 : double lifetime;
39 8523 : stream >> Z >> N >> decay.channel >> lifetime;
40 8523 : decay.rate = 1. / lifetime / c_light; // decay rate in [1/m]
41 : std::vector<double> gamma;
42 : double val;
43 34227 : while (stream >> val)
44 25704 : gamma.push_back(val);
45 21375 : for (int i = 0; i < gamma.size(); i += 2) {
46 12852 : decay.energy.push_back(gamma[i] * keV);
47 12852 : decay.intensity.push_back(gamma[i+1]);
48 : }
49 8523 : if (infile)
50 8523 : decayTable[Z * 31 + N].push_back(decay);
51 8541 : }
52 9 : infile.close();
53 18 : }
54 :
55 2 : void NuclearDecay::setHaveElectrons(bool b) {
56 2 : haveElectrons = b;
57 2 : }
58 :
59 1 : void NuclearDecay::setHavePhotons(bool b) {
60 1 : havePhotons = b;
61 1 : }
62 :
63 1 : void NuclearDecay::setHaveNeutrinos(bool b) {
64 1 : haveNeutrinos = b;
65 1 : }
66 :
67 0 : void NuclearDecay::setLimit(double l) {
68 0 : limit = l;
69 0 : }
70 :
71 32605 : void NuclearDecay::process(Candidate *candidate) const {
72 : // the loop should be processed at least once for limiting the next step
73 32605 : double step = candidate->getCurrentStep();
74 32605 : double z = candidate->getRedshift();
75 : do {
76 : // check if nucleus
77 32669 : int id = candidate->current.getId();
78 32669 : if (not (isNucleus(id)))
79 : return;
80 :
81 32668 : int A = massNumber(id);
82 32668 : int Z = chargeNumber(id);
83 32668 : int N = A - Z;
84 :
85 : // check if particle can decay
86 32668 : const std::vector<DecayMode> &decays = decayTable[Z * 31 + N];
87 32668 : if (decays.size() == 0)
88 : return;
89 :
90 : // find interaction mode with minimum random decay distance
91 340 : Random &random = Random::instance();
92 : double randDistance = std::numeric_limits<double>::max();
93 : int channel;
94 : double totalRate = 0;
95 :
96 707 : for (size_t i = 0; i < decays.size(); i++) {
97 367 : double rate = decays[i].rate;
98 367 : rate /= candidate->current.getLorentzFactor(); // relativistic time dilation
99 367 : rate /= (1 + z); // rate per light travel distance -> rate per comoving distance
100 367 : totalRate += rate;
101 367 : double d = -log(random.rand()) / rate;
102 367 : if (d > randDistance)
103 27 : continue;
104 : randDistance = d;
105 340 : channel = decays[i].channel;
106 : }
107 :
108 : // check if interaction doesn't happen
109 340 : if (step < randDistance) {
110 : // limit next step to a fraction of the mean free path
111 276 : candidate->limitNextStep(limit / totalRate);
112 276 : return;
113 : }
114 :
115 : // interact and repeat with remaining step
116 64 : performInteraction(candidate, channel);
117 64 : step -= randDistance;
118 64 : } while (step > 0);
119 : }
120 :
121 1023 : void NuclearDecay::performInteraction(Candidate *candidate, int channel) const {
122 : // interpret decay channel
123 : int nBetaMinus = digit(channel, 10000);
124 : int nBetaPlus = digit(channel, 1000);
125 : int nAlpha = digit(channel, 100);
126 : int nProton = digit(channel, 10);
127 : int nNeutron = digit(channel, 1);
128 :
129 : // perform decays
130 1023 : if (havePhotons)
131 11 : gammaEmission(candidate,channel);
132 1391 : for (size_t i = 0; i < nBetaMinus; i++)
133 368 : betaDecay(candidate, false);
134 1191 : for (size_t i = 0; i < nBetaPlus; i++)
135 168 : betaDecay(candidate, true);
136 1048 : for (size_t i = 0; i < nAlpha; i++)
137 25 : nucleonEmission(candidate, 4, 2);
138 1359 : for (size_t i = 0; i < nProton; i++)
139 336 : nucleonEmission(candidate, 1, 1);
140 1363 : for (size_t i = 0; i < nNeutron; i++)
141 340 : nucleonEmission(candidate, 1, 0);
142 1023 : }
143 :
144 11 : void NuclearDecay::gammaEmission(Candidate *candidate, int channel) const {
145 11 : int id = candidate->current.getId();
146 11 : int Z = chargeNumber(id);
147 11 : int N = massNumber(id) - Z;
148 :
149 : // get photon energies and emission probabilities for decay channel
150 11 : const std::vector<DecayMode> &decays = decayTable[Z * 31 + N];
151 : size_t idecay = decays.size();
152 21 : while (idecay-- != 0) {
153 21 : if (decays[idecay].channel == channel)
154 : break;
155 : }
156 :
157 : const std::vector<double> &energy = decays[idecay].energy;
158 : const std::vector<double> &intensity = decays[idecay].intensity;
159 :
160 : // check if photon emission available
161 11 : if (energy.size() == 0)
162 0 : return;
163 :
164 11 : Random &random = Random::instance();
165 11 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
166 :
167 26 : for (int i = 0; i < energy.size(); ++i) {
168 : // check if photon of specific energy is emitted
169 15 : if (random.rand() > intensity[i])
170 7 : continue;
171 : // create secondary photon; boost to lab frame
172 8 : double cosTheta = 2 * random.rand() - 1;
173 8 : double E = energy[i] * candidate->current.getLorentzFactor() * (1. - cosTheta);
174 8 : candidate->addSecondary(22, E, pos, 1., interactionTag);
175 : }
176 : }
177 :
178 536 : void NuclearDecay::betaDecay(Candidate *candidate, bool isBetaPlus) const {
179 536 : double gamma = candidate->current.getLorentzFactor();
180 536 : int id = candidate->current.getId();
181 536 : int A = massNumber(id);
182 536 : int Z = chargeNumber(id);
183 :
184 : // beta- decay
185 : int electronId = 11; // electron
186 : int neutrinoId = -12; // anti-electron neutrino
187 : int dZ = 1;
188 : // beta+ decay
189 536 : if (isBetaPlus) {
190 : electronId = -11; // positron
191 : neutrinoId = 12; // electron neutrino
192 : dZ = -1;
193 : }
194 :
195 : // update candidate, nuclear recoil negligible
196 : try
197 : {
198 536 : candidate->current.setId(nucleusId(A, Z + dZ));
199 : }
200 0 : catch (std::runtime_error &e)
201 : {
202 0 : KISS_LOG_ERROR<< "Something went wrong in the NuclearDecay\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();
203 0 : throw;
204 0 : }
205 :
206 536 : candidate->current.setLorentzFactor(gamma);
207 :
208 536 : if (not (haveElectrons or haveNeutrinos))
209 524 : return;
210 :
211 : // Q-value of the decay, subtract total energy of emitted photons
212 12 : double m1 = nuclearMass(A, Z);
213 12 : double m2 = nuclearMass(A, Z+dZ);
214 12 : double Q = (m1 - m2 - mass_electron) * c_squared;
215 :
216 : // generate cdf of electron energy, neglecting Coulomb correction
217 : // see Basdevant, Fundamentals in Nuclear Physics, eq. (4.92)
218 : // This leads to deviations from theoretical expectations at low
219 : // primary energies.
220 : std::vector<double> energies;
221 : std::vector<double> densities; // cdf(E), unnormalized
222 :
223 12 : energies.reserve(51);
224 12 : densities.reserve(51);
225 :
226 : double me = mass_electron * c_squared;
227 12 : double cdf = 0;
228 624 : for (int i = 0; i <= 50; i++) {
229 612 : double E = me + i / 50. * Q;
230 612 : cdf += E * sqrt(E * E - me * me) * pow(Q + me - E, 2);
231 612 : energies.push_back(E);
232 612 : densities.push_back(cdf);
233 : }
234 :
235 : // draw random electron energy and angle
236 : // assumption of ultra-relativistic particles
237 : // leads to deviations from theoretical predictions
238 : // is not problematic for usual CRPropa energies E>~TeV
239 12 : Random &random = Random::instance();
240 12 : double E = interpolate(random.rand() * cdf, densities, energies);
241 12 : double p = sqrt(E * E - me * me); // p*c
242 12 : double cosTheta = 2 * random.rand() - 1;
243 :
244 : // boost to lab frame
245 12 : double Ee = gamma * (E - p * cosTheta);
246 12 : double Enu = gamma * (Q + me - E) * (1 + cosTheta); // pnu*c ~ Enu
247 :
248 12 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
249 12 : if (haveElectrons)
250 12 : candidate->addSecondary(electronId, Ee, pos, 1., interactionTag);
251 12 : if (haveNeutrinos)
252 10 : candidate->addSecondary(neutrinoId, Enu, pos, 1., interactionTag);
253 : }
254 :
255 701 : void NuclearDecay::nucleonEmission(Candidate *candidate, int dA, int dZ) const {
256 701 : Random &random = Random::instance();
257 701 : int id = candidate->current.getId();
258 701 : int A = massNumber(id);
259 701 : int Z = chargeNumber(id);
260 701 : double EpA = candidate->current.getEnergy() / double(A);
261 :
262 : try
263 : {
264 701 : candidate->current.setId(nucleusId(A - dA, Z - dZ));
265 : }
266 0 : catch (std::runtime_error &e)
267 : {
268 0 : KISS_LOG_ERROR<< "Something went wrong in the NuclearDecay\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();
269 0 : throw;
270 0 : }
271 :
272 701 : candidate->current.setEnergy(EpA * (A - dA));
273 701 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(),candidate->current.getPosition());
274 :
275 : try
276 : {
277 701 : candidate->addSecondary(nucleusId(dA, dZ), EpA * dA, pos, 1., interactionTag);
278 : }
279 0 : catch (std::runtime_error &e)
280 : {
281 0 : KISS_LOG_ERROR<< "Something went wrong in the NuclearDecay\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();
282 0 : throw;
283 0 : }
284 :
285 701 : }
286 :
287 0 : double NuclearDecay::meanFreePath(int id, double gamma) {
288 0 : if (not (isNucleus(id)))
289 : return std::numeric_limits<double>::max();
290 :
291 0 : int A = massNumber(id);
292 0 : int Z = chargeNumber(id);
293 0 : int N = A - Z;
294 :
295 : // check if particle can decay
296 0 : const std::vector<DecayMode> &decays = decayTable[Z * 31 + N];
297 0 : if (decays.size() == 0)
298 : return std::numeric_limits<double>::max();
299 :
300 : double totalRate = 0;
301 :
302 0 : for (size_t i = 0; i < decays.size(); i++) {
303 0 : double rate = decays[i].rate;
304 0 : rate /= gamma;
305 0 : totalRate += rate;
306 : }
307 :
308 0 : return 1. / totalRate;
309 : }
310 :
311 1 : void NuclearDecay::setInteractionTag(std::string tag) {
312 1 : interactionTag = tag;
313 1 : }
314 :
315 2 : std::string NuclearDecay::getInteractionTag() const {
316 2 : return interactionTag;
317 : }
318 :
319 : } // namespace crpropa
|
