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