Line data Source code
1 : #include "crpropa/module/SynchrotronRadiation.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 2 : SynchrotronRadiation::SynchrotronRadiation(ref_ptr<MagneticField> field, bool havePhotons, double thinning, int nSamples, double limit) {
12 2 : setField(field);
13 2 : setBrms(0);
14 2 : initSpectrum();
15 2 : setHavePhotons(havePhotons);
16 2 : setLimit(limit);
17 2 : setSecondaryThreshold(1e6 * eV);
18 2 : setMaximumSamples(nSamples);
19 2 : setThinning(thinning);
20 2 : }
21 :
22 6 : SynchrotronRadiation::SynchrotronRadiation(double Brms, bool havePhotons, double thinning, int nSamples, double limit) {
23 6 : setBrms(Brms);
24 6 : initSpectrum();
25 6 : setHavePhotons(havePhotons);
26 6 : setLimit(limit);
27 6 : setSecondaryThreshold(1e6 * eV);
28 6 : setMaximumSamples(nSamples);
29 6 : setThinning(thinning);
30 6 : }
31 :
32 3 : void SynchrotronRadiation::setField(ref_ptr<MagneticField> f) {
33 3 : this->field = f;
34 3 : }
35 :
36 3 : ref_ptr<MagneticField> SynchrotronRadiation::getField() {
37 3 : return field;
38 : }
39 :
40 9 : void SynchrotronRadiation::setBrms(double Brms) {
41 9 : this->Brms = Brms;
42 9 : }
43 :
44 5 : double SynchrotronRadiation::getBrms() {
45 5 : return Brms;
46 : }
47 :
48 9 : void SynchrotronRadiation::setHavePhotons(bool havePhotons) {
49 9 : this->havePhotons = havePhotons;
50 9 : }
51 :
52 5 : bool SynchrotronRadiation::getHavePhotons() {
53 5 : return havePhotons;
54 : }
55 :
56 9 : void SynchrotronRadiation::setThinning(double thinning) {
57 9 : this->thinning = thinning;
58 9 : }
59 :
60 5 : double SynchrotronRadiation::getThinning() {
61 5 : return thinning;
62 : }
63 :
64 9 : void SynchrotronRadiation::setLimit(double limit) {
65 9 : this->limit = limit;
66 9 : }
67 :
68 5 : double SynchrotronRadiation::getLimit() {
69 5 : return limit;
70 : }
71 :
72 10 : void SynchrotronRadiation::setMaximumSamples(int nmax) {
73 10 : maximumSamples = nmax;
74 10 : }
75 :
76 5 : int SynchrotronRadiation::getMaximumSamples() {
77 5 : return maximumSamples;
78 : }
79 :
80 10 : void SynchrotronRadiation::setSecondaryThreshold(double threshold) {
81 10 : secondaryThreshold = threshold;
82 10 : }
83 :
84 5 : double SynchrotronRadiation::getSecondaryThreshold() const {
85 5 : return secondaryThreshold;
86 : }
87 :
88 8 : void SynchrotronRadiation::initSpectrum() {
89 16 : std::string filename = getDataPath("Synchrotron/spectrum.txt");
90 8 : std::ifstream infile(filename.c_str());
91 :
92 8 : if (!infile.good())
93 0 : throw std::runtime_error("SynchrotronRadiation: could not open file " + filename);
94 :
95 : // clear previously loaded interaction rates
96 : tabx.clear();
97 : tabCDF.clear();
98 :
99 11256 : while (infile.good()) {
100 11248 : if (infile.peek() != '#') {
101 : double a, b;
102 : infile >> a >> b;
103 11216 : if (infile) {
104 11208 : tabx.push_back(pow(10, a));
105 11208 : tabCDF.push_back(b);
106 : }
107 : }
108 11248 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
109 : }
110 8 : infile.close();
111 16 : }
112 :
113 6 : void SynchrotronRadiation::process(Candidate *candidate) const {
114 6 : double charge = fabs(candidate->current.getCharge());
115 6 : if (charge == 0)
116 4 : return; // only charged particles
117 :
118 : // calculate gyroradius, evaluated at the current position
119 6 : double z = candidate->getRedshift();
120 : double B;
121 6 : if (field.valid()) {
122 0 : Vector3d Bvec = field->getField(candidate->current.getPosition(), z);
123 0 : B = Bvec.cross(candidate->current.getDirection()).getR();
124 : } else {
125 6 : B = sqrt(2. / 3) * Brms; // average perpendicular field component
126 : }
127 6 : B *= pow(1 + z, 2); // cosmological scaling
128 6 : double Rg = candidate->current.getMomentum().getR() / charge / B;
129 :
130 : // calculate energy loss
131 6 : double lf = candidate->current.getLorentzFactor();
132 6 : double dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
133 6 : double step = candidate->getCurrentStep() / (1 + z); // step size in local frame
134 6 : double dE = step * dEdx;
135 :
136 : // apply energy loss and limit next step
137 6 : double E = candidate->current.getEnergy();
138 6 : candidate->current.setEnergy(E - dE);
139 6 : candidate->limitNextStep(limit * E / dEdx);
140 :
141 : // optionally add secondary photons
142 6 : if (not(havePhotons))
143 : return;
144 :
145 : // check if photons with energies > 14 * Ecrit are possible
146 2 : double Ecrit = 3. / 4 * h_planck / M_PI * c_light * pow(lf, 3) / Rg;
147 2 : if (14 * Ecrit < secondaryThreshold)
148 : return;
149 :
150 : // draw photons up to the total energy loss
151 : // if maximumSamples is reached before that, compensate the total energy afterwards
152 2 : Random &random = Random::instance();
153 : double dE0 = dE;
154 : std::vector<double> energies;
155 : int counter = 0;
156 3626394 : while (dE > 0) {
157 : // draw random value between 0 and maximum of corresponding cdf
158 : // choose bin of s where cdf(x) = cdf_rand -> x_rand
159 3626394 : size_t i = random.randBin(tabCDF); // draw random bin (upper bin boundary returned)
160 3626394 : double binWidth = (tabx[i] - tabx[i-1]);
161 3626394 : double x = tabx[i-1] + random.rand() * binWidth; // draw random x uniformly distributed in bin
162 3626394 : double Ephoton = x * Ecrit;
163 :
164 : // if the remaining energy is not sufficient check for random accepting
165 3626394 : if (Ephoton > dE) {
166 1 : if (random.rand() > (dE / Ephoton))
167 : break; // not accepted
168 : }
169 :
170 : // only activate the "per-step" sampling if maximumSamples is explicitly set.
171 3626393 : if (maximumSamples > 0) {
172 1001 : if (counter >= maximumSamples)
173 : break;
174 : }
175 :
176 : // store energies in array
177 3626392 : energies.push_back(Ephoton);
178 :
179 : // energy loss
180 3626392 : dE -= Ephoton;
181 :
182 : // counter for sampling break condition;
183 3626392 : counter++;
184 : }
185 :
186 : // while loop before gave total energy which is just a fraction of the required
187 : double w1 = 1;
188 2 : if (maximumSamples > 0 && dE > 0)
189 1 : w1 = 1. / (1. - dE / dE0);
190 :
191 : // loop over sampled photons and attribute weights accordingly
192 3626394 : for (int i = 0; i < energies.size(); i++) {
193 3626392 : double Ephoton = energies[i];
194 3626392 : double f = Ephoton / (E - dE0);
195 3626392 : double w = w1 / pow(f, thinning);
196 :
197 : // thinning procedure: accepts only a few random secondaries
198 3626392 : if (random.rand() < pow(f, thinning)) {
199 3626392 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
200 3626392 : if (Ephoton > secondaryThreshold) // create only photons with energies above threshold
201 3314271 : candidate->addSecondary(22, Ephoton, pos, w, interactionTag);
202 : }
203 : }
204 2 : }
205 :
206 0 : std::string SynchrotronRadiation::getDescription() const {
207 0 : std::stringstream s;
208 0 : s << "Synchrotron radiation";
209 0 : if (field.valid())
210 0 : s << " for specified magnetic field";
211 : else
212 0 : s << " for Brms = " << Brms / nG << " nG";
213 0 : if (havePhotons)
214 0 : s << ", synchrotron photons E > " << secondaryThreshold / eV << " eV";
215 : else
216 0 : s << ", no synchrotron photons";
217 0 : if (maximumSamples > 0)
218 0 : s << "maximum number of photon samples: " << maximumSamples;
219 0 : if (thinning > 0)
220 0 : s << "thinning parameter: " << thinning;
221 0 : return s.str();
222 0 : }
223 :
224 1 : void SynchrotronRadiation::setInteractionTag(std::string tag) {
225 1 : interactionTag = tag;
226 1 : }
227 :
228 2 : std::string SynchrotronRadiation::getInteractionTag() const {
229 2 : return interactionTag;
230 : }
231 :
232 : } // namespace crpropa
|
