Line data Source code
1 : #include "crpropa/module/ElectronPairProduction.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 <stdexcept>
10 :
11 : namespace crpropa {
12 :
13 10 : ElectronPairProduction::ElectronPairProduction(ref_ptr<PhotonField> photonField,
14 10 : bool haveElectrons, double limit) {
15 10 : this->haveElectrons = haveElectrons;
16 10 : this->limit = limit;
17 10 : setPhotonField(photonField);
18 10 : }
19 :
20 19 : void ElectronPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
21 19 : this->photonField = photonField;
22 19 : std::string fname = photonField->getFieldName();
23 19 : setDescription("ElectronPairProduction: " + fname);
24 38 : initRate(getDataPath("ElectronPairProduction/lossrate_" + fname + ".txt"));
25 19 : if (haveElectrons) { // Load secondary spectra only if electrons should be produced
26 0 : initSpectrum(getDataPath("ElectronPairProduction/spectrum_" + fname.substr(0,3) + ".txt"));
27 : }
28 19 : }
29 :
30 1 : void ElectronPairProduction::setHaveElectrons(bool haveElectrons) {
31 1 : this->haveElectrons = haveElectrons;
32 1 : if (haveElectrons) { // Load secondary spectra in case haveElectrons was changed to true
33 1 : std::string fname = photonField->getFieldName();
34 3 : initSpectrum(getDataPath("ElectronPairProduction/spectrum_" + fname.substr(0,3) + ".txt"));
35 : }
36 1 : }
37 :
38 0 : void ElectronPairProduction::setLimit(double limit) {
39 0 : this->limit = limit;
40 0 : }
41 :
42 19 : void ElectronPairProduction::initRate(std::string filename) {
43 19 : std::ifstream infile(filename.c_str());
44 :
45 19 : if (!infile.good())
46 0 : throw std::runtime_error("ElectronPairProduction: could not open file " + filename);
47 :
48 : // clear previously loaded interaction rates
49 : tabLorentzFactor.clear();
50 : tabLossRate.clear();
51 :
52 2853 : while (infile.good()) {
53 2834 : if (infile.peek() != '#') {
54 : double a, b;
55 : infile >> a >> b;
56 2777 : if (infile) {
57 2758 : tabLorentzFactor.push_back(pow(10, a));
58 2758 : tabLossRate.push_back(b / Mpc);
59 : }
60 : }
61 2834 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
62 : }
63 19 : infile.close();
64 19 : }
65 :
66 1 : void ElectronPairProduction::initSpectrum(std::string filename) {
67 1 : std::ifstream infile(filename.c_str());
68 1 : if (!infile.good())
69 0 : throw std::runtime_error("ElectronPairProduction: could not open file " + filename);
70 :
71 : double dNdE;
72 1 : tabSpectrum.resize(70);
73 71 : for (size_t i = 0; i < 70; i++) {
74 70 : tabSpectrum[i].resize(170);
75 11970 : for (size_t j = 0; j < 170; j++) {
76 : infile >> dNdE;
77 11900 : tabSpectrum[i][j] = dNdE * pow(10, (7 + 0.1 * j)); // read electron distribution pdf(Ee) ~ dN/dEe * Ee
78 : }
79 11900 : for (size_t j = 1; j < 170; j++) {
80 11830 : tabSpectrum[i][j] += tabSpectrum[i][j - 1]; // cdf(Ee), unnormalized
81 : }
82 : }
83 1 : infile.close();
84 1 : }
85 :
86 1030 : double ElectronPairProduction::lossLength(int id, double lf, double z) const {
87 1030 : double Z = chargeNumber(id);
88 1030 : if (Z == 0)
89 : return std::numeric_limits<double>::max(); // no pair production on uncharged particles
90 :
91 1018 : lf *= (1 + z);
92 1018 : if (lf < tabLorentzFactor.front())
93 : return std::numeric_limits<double>::max(); // below energy threshold
94 :
95 : double rate;
96 995 : if (lf < tabLorentzFactor.back())
97 995 : rate = interpolate(lf, tabLorentzFactor, tabLossRate); // interpolation
98 : else
99 0 : rate = tabLossRate.back() * pow(lf / tabLorentzFactor.back(), -0.6); // extrapolation
100 :
101 995 : double A = nuclearMass(id) / mass_proton; // more accurate than massNumber(Id)
102 995 : rate *= Z * Z / A * pow_integer<3>(1 + z) * photonField->getRedshiftScaling(z);
103 995 : return 1. / rate;
104 : }
105 :
106 1031 : void ElectronPairProduction::process(Candidate *c) const {
107 1031 : int id = c->current.getId();
108 1031 : if (not (isNucleus(id)))
109 : return; // only nuclei
110 :
111 1030 : double lf = c->current.getLorentzFactor();
112 1030 : double z = c->getRedshift();
113 1030 : double losslen = lossLength(id, lf, z); // energy loss length
114 1030 : if (losslen >= std::numeric_limits<double>::max())
115 : return;
116 :
117 995 : double step = c->getCurrentStep() / (1 + z); // step size in local frame
118 995 : double loss = step / losslen; // relative energy loss
119 :
120 995 : if (haveElectrons) {
121 1 : double dE = c->current.getEnergy() * loss; // energy loss
122 1 : int i = round((log10(lf) - 6.05) * 10); // find closest cdf(Ee|log10(gamma))
123 1 : i = std::min(std::max(i, 0), 69);
124 1 : Random &random = Random::instance();
125 :
126 : // draw pairs as long as their energy is smaller than the pair production energy loss
127 6173 : while (dE > 0) {
128 6172 : size_t j = random.randBin(tabSpectrum[i]);
129 6172 : double Ee = pow(10, 6.95 + (j + random.rand()) * 0.1) * eV;
130 6172 : double Epair = 2 * Ee; // NOTE: electron and positron in general don't have same lab frame energy, but averaged over many draws the result is consistent
131 : // if the remaining energy is not sufficient check for random accepting
132 6172 : if (Epair > dE)
133 1 : if (random.rand() > (dE / Epair))
134 : break; // not accepted
135 :
136 : // create pair and repeat with remaining energy
137 6172 : dE -= Epair;
138 6172 : Vector3d pos = random.randomInterpolatedPosition(c->previous.getPosition(), c->current.getPosition());
139 6172 : c->addSecondary( 11, Ee, pos, 1., interactionTag);
140 6172 : c->addSecondary(-11, Ee, pos, 1., interactionTag);
141 : }
142 : }
143 :
144 995 : c->current.setLorentzFactor(lf * (1 - loss));
145 995 : c->limitNextStep(limit * losslen);
146 : }
147 :
148 1 : void ElectronPairProduction::setInteractionTag(std::string tag) {
149 1 : interactionTag = tag;
150 1 : }
151 :
152 2 : std::string ElectronPairProduction::getInteractionTag() const {
153 2 : return interactionTag;
154 : }
155 :
156 :
157 : } // namespace crpropa
|
