Line data Source code
1 : #include "crpropa/module/EMTripletPairProduction.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 :
5 : #include <fstream>
6 : #include <limits>
7 : #include <stdexcept>
8 :
9 : namespace crpropa {
10 :
11 : static const double mec2 = mass_electron * c_squared;
12 :
13 3 : EMTripletPairProduction::EMTripletPairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit) {
14 3 : setPhotonField(photonField);
15 3 : setHaveElectrons(haveElectrons);
16 3 : setLimit(limit);
17 3 : setThinning(thinning);
18 3 : }
19 :
20 15 : void EMTripletPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
21 15 : this->photonField = photonField;
22 15 : std::string fname = photonField->getFieldName();
23 15 : setDescription("EMTripletPairProduction: " + fname);
24 30 : initRate(getDataPath("EMTripletPairProduction/rate_" + fname + ".txt"));
25 45 : initCumulativeRate(getDataPath("EMTripletPairProduction/cdf_" + fname + ".txt"));
26 15 : }
27 :
28 4 : void EMTripletPairProduction::setHaveElectrons(bool haveElectrons) {
29 4 : this->haveElectrons = haveElectrons;
30 4 : }
31 :
32 3 : void EMTripletPairProduction::setLimit(double limit) {
33 3 : this->limit = limit;
34 3 : }
35 :
36 3 : void EMTripletPairProduction::setThinning(double thinning) {
37 3 : this->thinning = thinning;
38 3 : }
39 :
40 15 : void EMTripletPairProduction::initRate(std::string filename) {
41 15 : std::ifstream infile(filename.c_str());
42 :
43 15 : if (!infile.good())
44 0 : throw std::runtime_error("EMTripletPairProduction: could not open file " + filename);
45 :
46 : // clear previously loaded interaction rates
47 : tabEnergy.clear();
48 : tabRate.clear();
49 :
50 3340 : while (infile.good()) {
51 3325 : if (infile.peek() != '#') {
52 : double a, b;
53 : infile >> a >> b;
54 3265 : if (infile) {
55 3250 : tabEnergy.push_back(pow(10, a) * eV);
56 3250 : tabRate.push_back(b / Mpc);
57 : }
58 : }
59 3325 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
60 : }
61 15 : infile.close();
62 15 : }
63 :
64 15 : void EMTripletPairProduction::initCumulativeRate(std::string filename) {
65 15 : std::ifstream infile(filename.c_str());
66 :
67 15 : if (!infile.good())
68 : throw std::runtime_error(
69 0 : "EMTripletPairProduction: could not open file " + filename);
70 :
71 : // clear previously loaded tables
72 : tabE.clear();
73 : tabs.clear();
74 : tabCDF.clear();
75 :
76 : // skip header
77 75 : while (infile.peek() == '#')
78 60 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
79 :
80 : // read s values in first line
81 : double a;
82 : infile >> a; // skip first value
83 1590 : while (infile.good() and (infile.peek() != '\n')) {
84 : infile >> a;
85 1575 : tabs.push_back(pow(10, a) * eV * eV);
86 : }
87 :
88 : // read all following lines: E, cdf values
89 3265 : while (infile.good()) {
90 : infile >> a;
91 3265 : if (!infile)
92 : break; // end of file
93 3250 : tabE.push_back(pow(10, a) * eV);
94 : std::vector<double> cdf;
95 344500 : for (int i = 0; i < tabs.size(); i++) {
96 : infile >> a;
97 341250 : cdf.push_back(a / Mpc);
98 : }
99 3250 : tabCDF.push_back(cdf);
100 3250 : }
101 15 : infile.close();
102 15 : }
103 :
104 1 : void EMTripletPairProduction::performInteraction(Candidate *candidate) const {
105 1 : int id = candidate->current.getId();
106 1 : if (abs(id) != 11)
107 : return;
108 :
109 : // scale the particle energy instead of background photons
110 1 : double z = candidate->getRedshift();
111 1 : double E = candidate->current.getEnergy() * (1 + z);
112 :
113 1 : if (E < tabE.front() or E > tabE.back())
114 : return;
115 :
116 : // sample the value of eps
117 1 : Random &random = Random::instance();
118 1 : size_t i = closestIndex(E, tabE);
119 1 : size_t j = random.randBin(tabCDF[i]);
120 1 : double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1);
121 1 : double eps = s_kin / 4. / E; // random background photon energy
122 :
123 : // Use approximation from A. Mastichiadis et al., Astroph. Journ. 300:178-189 (1986), eq. 30.
124 : // This approx is valid only for alpha >=100 where alpha = p0*eps*costheta - E0*eps
125 : // For our purposes, me << E0 --> p0~E0 --> alpha = E0*eps*(costheta - 1) >= 100
126 1 : double Epp = 5.7e-1 * pow(eps / mec2, -0.56) * pow(E / mec2, 0.44) * mec2;
127 :
128 1 : double f = Epp / E;
129 :
130 1 : if (haveElectrons) {
131 1 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
132 1 : if (random.rand() < pow(1 - f, thinning)) {
133 1 : double w = 1. / pow(1 - f, thinning);
134 1 : candidate->addSecondary(11, Epp / (1 + z), pos, w, interactionTag);
135 : }
136 1 : if (random.rand() < pow(f, thinning)) {
137 1 : double w = 1. / pow(f, thinning);
138 1 : candidate->addSecondary(-11, Epp / (1 + z), pos, w, interactionTag);
139 : }
140 : }
141 : // Update the primary particle energy.
142 : // This is done after adding the secondaries to correctly set the secondaries parent
143 1 : candidate->current.setEnergy((E - 2 * Epp) / (1. + z));
144 : }
145 :
146 1 : void EMTripletPairProduction::process(Candidate *candidate) const {
147 : // check if electron / positron
148 1 : int id = candidate->current.getId();
149 1 : if (abs(id) != 11)
150 : return;
151 :
152 : // scale the particle energy instead of background photons
153 1 : double z = candidate->getRedshift();
154 1 : double E = (1 + z) * candidate->current.getEnergy();
155 :
156 : // check if in tabulated energy range
157 1 : if ((E < tabEnergy.front()) or (E > tabEnergy.back()))
158 : return;
159 :
160 : // cosmological scaling of interaction distance (comoving)
161 1 : double scaling = pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
162 1 : double rate = scaling * interpolate(E, tabEnergy, tabRate);
163 :
164 : // run this loop at least once to limit the step size
165 1 : double step = candidate->getCurrentStep();
166 1 : Random &random = Random::instance();
167 : do {
168 1 : double randDistance = -log(random.rand()) / rate;
169 : // check for interaction; if it doesn't occur, limit next step
170 1 : if (step < randDistance) {
171 1 : candidate->limitNextStep(limit / rate);
172 1 : return;
173 : }
174 0 : performInteraction(candidate);
175 0 : step -= randDistance;
176 0 : } while (step > 0.);
177 : }
178 :
179 1 : void EMTripletPairProduction::setInteractionTag(std::string tag) {
180 1 : interactionTag = tag;
181 1 : }
182 :
183 2 : std::string EMTripletPairProduction::getInteractionTag() const {
184 2 : return interactionTag;
185 : }
186 :
187 : } // namespace crpropa
|
