Line data Source code
1 : #include "crpropa/Candidate.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/ParticleID.h"
4 : #include "crpropa/PhotonBackground.h"
5 : #include "crpropa/module/ElectronPairProduction.h"
6 : #include "crpropa/module/NuclearDecay.h"
7 : #include "crpropa/module/PhotoDisintegration.h"
8 : #include "crpropa/module/ElasticScattering.h"
9 : #include "crpropa/module/PhotoPionProduction.h"
10 : #include "crpropa/module/Redshift.h"
11 : #include "crpropa/module/EMPairProduction.h"
12 : #include "crpropa/module/EMDoublePairProduction.h"
13 : #include "crpropa/module/EMTripletPairProduction.h"
14 : #include "crpropa/module/EMInverseComptonScattering.h"
15 : #include "crpropa/module/SynchrotronRadiation.h"
16 : #include "gtest/gtest.h"
17 : #include "kiss/path.h"
18 :
19 : #include <fstream>
20 :
21 : namespace crpropa {
22 :
23 : // ElectronPairProduction -----------------------------------------------------
24 1 : TEST(ElectronPairProduction, allBackgrounds) {
25 : // Test if interaction data files are loaded.
26 1 : ref_ptr<PhotonField> cmb = new CMB();
27 1 : ElectronPairProduction epp(cmb);
28 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
29 1 : epp.setPhotonField(irb);
30 1 : irb = new IRB_Stecker05();
31 1 : epp.setPhotonField(irb);
32 1 : irb = new IRB_Franceschini08();
33 1 : epp.setPhotonField(irb);
34 1 : irb = new IRB_Finke10();
35 1 : epp.setPhotonField(irb);
36 1 : irb = new IRB_Dominguez11();
37 1 : epp.setPhotonField(irb);
38 1 : irb = new IRB_Gilmore12();
39 1 : epp.setPhotonField(irb);
40 1 : irb = new IRB_Stecker16_upper();
41 1 : epp.setPhotonField(irb);
42 1 : irb = new IRB_Stecker16_lower();
43 1 : epp.setPhotonField(irb);
44 1 : irb = new IRB_Finke22();
45 2 : epp.setPhotonField(irb);
46 2 : }
47 :
48 1 : TEST(ElectronPairProduction, energyDecreasing) {
49 : // Test if energy loss occurs for protons with energies from 1e15 - 1e23 eV.
50 1 : Candidate c;
51 1 : c.setCurrentStep(2 * Mpc);
52 1 : c.current.setId(nucleusId(1, 1)); // proton
53 :
54 1 : ref_ptr<PhotonField> cmb = new CMB();
55 1 : ElectronPairProduction epp1(cmb);
56 81 : for (int i = 0; i < 80; i++) {
57 80 : double E = pow(10, 15 + i * 0.1) * eV;
58 80 : c.current.setEnergy(E);
59 80 : epp1.process(&c);
60 80 : EXPECT_LE(c.current.getEnergy(), E);
61 : }
62 :
63 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
64 2 : ElectronPairProduction epp2(irb);
65 81 : for (int i = 0; i < 80; i++) {
66 80 : double E = pow(10, 15 + i * 0.1) * eV;
67 80 : c.current.setEnergy(E);
68 80 : epp2.process(&c);
69 80 : EXPECT_LE(c.current.getEnergy(), E);
70 : }
71 2 : }
72 :
73 1 : TEST(ElectronPairProduction, belowEnergyTreshold) {
74 : // Test if nothing happens below 1e15 eV.
75 1 : ref_ptr<PhotonField> cmb = new CMB();
76 1 : ElectronPairProduction epp(cmb);
77 1 : Candidate c(nucleusId(1, 1), 1E14 * eV);
78 1 : epp.process(&c);
79 1 : EXPECT_DOUBLE_EQ(1E14 * eV, c.current.getEnergy());
80 2 : }
81 :
82 1 : TEST(ElectronPairProduction, thisIsNotNucleonic) {
83 : // Test if non-nuclei are skipped.
84 1 : ref_ptr<PhotonField> cmb = new CMB();
85 1 : ElectronPairProduction epp(cmb);
86 1 : Candidate c(11, 1E20 * eV); // electron
87 1 : epp.process(&c);
88 1 : EXPECT_DOUBLE_EQ(1E20 * eV, c.current.getEnergy());
89 2 : }
90 :
91 2 : TEST(ElectronPairProduction, valuesCMB) {
92 : // Test if energy loss corresponds to the data table.
93 : std::vector<double> x;
94 : std::vector<double> y;
95 1 : std::ifstream infile(getDataPath("pair_CMB.txt").c_str());
96 1 : while (infile.good()) {
97 0 : if (infile.peek() != '#') {
98 : double a, b;
99 : infile >> a >> b;
100 0 : if (infile) {
101 0 : x.push_back(a * eV);
102 0 : y.push_back(b * eV / Mpc);
103 : }
104 : }
105 0 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
106 : }
107 1 : infile.close();
108 :
109 1 : Candidate c;
110 1 : c.setCurrentStep(1 * Mpc);
111 1 : c.current.setId(nucleusId(1, 1)); // proton
112 1 : ref_ptr<PhotonField> cmb = new CMB();
113 :
114 1 : ElectronPairProduction epp(cmb);
115 1 : for (int i = 0; i < x.size(); i++) {
116 0 : c.current.setEnergy(x[i]);
117 0 : epp.process(&c);
118 0 : double dE = x[i] - c.current.getEnergy();
119 0 : double dE_table = y[i] * 1 * Mpc;
120 0 : EXPECT_NEAR(dE_table, dE, 1e-12);
121 : }
122 2 : }
123 :
124 1 : TEST(ElectronPairProduction, interactionTag) {
125 :
126 1 : ref_ptr<PhotonField> cmb = new CMB();
127 1 : ElectronPairProduction epp(cmb);
128 :
129 : // test the default interaction tag
130 1 : EXPECT_TRUE(epp.getInteractionTag() == "EPP");
131 :
132 : // test changing the interaction tag
133 1 : epp.setInteractionTag("myTag");
134 1 : EXPECT_TRUE(epp.getInteractionTag() == "myTag");
135 :
136 : // test the tag of produced secondaries
137 2 : Candidate c;
138 1 : c.setCurrentStep(1 * Gpc);
139 1 : c.current.setId(nucleusId(1,1));
140 1 : c.current.setEnergy(100 * EeV);
141 1 : epp.setHaveElectrons(true);
142 1 : epp.process(&c);
143 :
144 1 : std::string secondaryTag = c.secondaries[0] -> getTagOrigin();
145 1 : EXPECT_TRUE(secondaryTag == "myTag");
146 2 : }
147 :
148 2 : TEST(ElectronPairProduction, valuesIRB) {
149 : // Test if energy loss corresponds to the data table.
150 : std::vector<double> x;
151 : std::vector<double> y;
152 1 : std::ifstream infile(getDataPath("pairIRB.txt").c_str());
153 1 : while (infile.good()) {
154 0 : if (infile.peek() != '#') {
155 : double a, b;
156 : infile >> a >> b;
157 0 : if (infile) {
158 0 : x.push_back(a * eV);
159 0 : y.push_back(b * eV / Mpc);
160 : }
161 : }
162 0 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
163 : }
164 1 : infile.close();
165 :
166 1 : Candidate c;
167 1 : c.setCurrentStep(1 * Mpc);
168 1 : c.current.setId(nucleusId(1, 1)); // proton
169 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
170 :
171 1 : ElectronPairProduction epp(irb);
172 1 : for (int i = 0; i < x.size(); i++) {
173 0 : c.current.setEnergy(x[i]);
174 0 : epp.process(&c);
175 0 : double dE = x[i] - c.current.getEnergy();
176 0 : double dE_table = y[i] * 1 * Mpc;
177 0 : EXPECT_NEAR(dE, dE_table, 1e-12);
178 : }
179 2 : }
180 :
181 : // NuclearDecay ---------------------------------------------------------------
182 1 : TEST(NuclearDecay, scandium44) {
183 : // Test beta+ decay of 44Sc to 44Ca.
184 : // This test can stochastically fail.
185 1 : NuclearDecay d(true, true);
186 1 : Candidate c(nucleusId(44, 21), 1E18 * eV);
187 1 : c.setCurrentStep(100 * Mpc);
188 1 : double gamma = c.current.getLorentzFactor();
189 1 : d.process(&c);
190 :
191 : // expected decay product: 44Ca
192 1 : EXPECT_EQ(nucleusId(44, 20), c.current.getId());
193 :
194 : // expect Lorentz factor to be conserved
195 1 : EXPECT_DOUBLE_EQ(gamma, c.current.getLorentzFactor());
196 :
197 : // expect at least two secondaries: positron + electron neutrino
198 1 : EXPECT_GE(c.secondaries.size(), 2);
199 1 : }
200 :
201 1 : TEST(NuclearDecay, lithium4) {
202 : // Test proton dripping of Li-4 to He-3
203 : // This test can stochastically fail
204 1 : NuclearDecay d;
205 1 : Candidate c(nucleusId(4, 3), 4 * EeV);
206 1 : c.setCurrentStep(100 * Mpc);
207 1 : d.process(&c);
208 :
209 : // expected decay product: He-3
210 1 : EXPECT_EQ(nucleusId(3, 2), c.current.getId());
211 :
212 : // expected secondary: proton
213 1 : EXPECT_EQ(1, c.secondaries.size());
214 2 : Candidate c1 = *c.secondaries[0];
215 1 : EXPECT_EQ(nucleusId(1, 1), c1.current.getId());
216 1 : EXPECT_EQ(1 * EeV, c1.current.getEnergy());
217 1 : }
218 :
219 1 : TEST(NuclearDecay, helium5) {
220 : // Test neutron dripping of He-5 to He-4.
221 : // This test can stochastically fail.
222 1 : NuclearDecay d;
223 1 : Candidate c(nucleusId(5, 2), 5 * EeV);
224 1 : c.setCurrentStep(100 * Mpc);
225 1 : d.process(&c);
226 :
227 : // expected primary: He-4
228 1 : EXPECT_EQ(nucleusId(4, 2), c.current.getId());
229 1 : EXPECT_EQ(4, c.current.getEnergy() / EeV);
230 :
231 : // expected secondary: neutron
232 2 : Candidate c2 = *c.secondaries[0];
233 1 : EXPECT_EQ(nucleusId(1, 0), c2.current.getId());
234 1 : EXPECT_EQ(1, c2.current.getEnergy() / EeV);
235 1 : }
236 :
237 1 : TEST(NuclearDecay, limitNextStep) {
238 : // Test if next step is limited in case of a neutron.
239 1 : NuclearDecay decay;
240 1 : Candidate c(nucleusId(1, 0), 10 * EeV);
241 1 : c.setNextStep(std::numeric_limits<double>::max());
242 1 : decay.process(&c);
243 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
244 1 : }
245 :
246 1 : TEST(NuclearDecay, allChannelsWorking) {
247 : // Test if all nuclear decays are working.
248 1 : NuclearDecay d;
249 1 : Candidate c;
250 :
251 1 : std::ifstream infile(getDataPath("nuclear_decay.txt").c_str());
252 11556 : while (infile.good()) {
253 11555 : if (infile.peek() != '#') {
254 : int Z, N, channel, foo;
255 11553 : infile >> Z >> N >> channel >> foo;
256 11553 : c.current.setId(nucleusId(Z + N, Z));
257 11553 : c.current.setEnergy(80 * EeV);
258 11553 : d.performInteraction(&c, channel);
259 : }
260 11555 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
261 : }
262 1 : infile.close();
263 1 : }
264 :
265 1 : TEST(NuclearDecay, secondaries) {
266 : // Test if all types of secondaries are produced.
267 1 : NuclearDecay d;
268 1 : d.setHaveElectrons(true);
269 1 : d.setHaveNeutrinos(true);
270 1 : d.setHavePhotons(true);
271 1 : Candidate c;
272 :
273 : // He-8 --> Li-8 + e- + neutrino
274 : // additional photon emitted with 84% probability
275 : // --> expect at least 1 photon out of 10 decays
276 11 : for (int i = 0; i < 10; ++i) {
277 10 : c.current.setId(nucleusId(8, 2));
278 10 : c.current.setEnergy(5 * EeV);
279 10 : d.performInteraction(&c, 10000);
280 : }
281 :
282 : // count number of secondaries
283 1 : size_t nElectrons = 0;
284 1 : size_t nNeutrinos = 0;
285 1 : size_t nPhotons = 0;
286 :
287 27 : for(size_t i = 0; i < c.secondaries.size(); ++i) {
288 26 : int id = (*c.secondaries[i]).current.getId();
289 26 : if (id == 22) nPhotons++;
290 26 : if (id == 11) nElectrons++;
291 26 : if (id == -12) nNeutrinos++;
292 : }
293 :
294 1 : EXPECT_EQ(nElectrons, 10);
295 1 : EXPECT_EQ(nNeutrinos, 10);
296 1 : EXPECT_GE(nPhotons, 1);
297 1 : }
298 :
299 1 : TEST(NuclearDecay, thisIsNotNucleonic) {
300 : // Test if nothing happens to an electron
301 1 : NuclearDecay decay;
302 1 : Candidate c(11, 10 * EeV);
303 1 : c.setNextStep(std::numeric_limits<double>::max());
304 1 : decay.process(&c);
305 1 : EXPECT_EQ(11, c.current.getId());
306 1 : EXPECT_EQ(10 * EeV, c.current.getEnergy());
307 1 : }
308 :
309 1 : TEST(NuclearDecay, interactionTag) {
310 : // test default interaction tag
311 1 : NuclearDecay decay;
312 1 : EXPECT_TRUE(decay.getInteractionTag() == "ND");
313 :
314 : // test secondary tag
315 1 : decay.setHaveElectrons(true);
316 2 : Candidate c(nucleusId(8,2), 5 * EeV);
317 1 : decay.performInteraction(&c, 10000);
318 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "ND");
319 :
320 : // test custom tags
321 1 : decay.setInteractionTag("myTag");
322 1 : EXPECT_TRUE(decay.getInteractionTag() == "myTag");
323 1 : }
324 :
325 1 : TEST(NuclearDecay, superheavy_stableHeavyNucleusIsNoOp) {
326 : // Zn-68 (Z=30, N=38) is a stable nucleus and has no entry in
327 : // nuclear_decay.txt even in the extended (Z<=82) table.
328 : // process() must hit the empty-vector guard and return without
329 : // touching the candidate.
330 1 : NuclearDecay d;
331 1 : Candidate c(nucleusId(68, 30), 80 * EeV); // Zn-68, Z=30
332 1 : c.setCurrentStep(100 * Mpc);
333 1 : int id_before = c.current.getId();
334 1 : double E_before = c.current.getEnergy();
335 1 : d.process(&c);
336 1 : EXPECT_EQ(id_before, c.current.getId());
337 1 : EXPECT_DOUBLE_EQ(E_before, c.current.getEnergy());
338 1 : EXPECT_TRUE(c.isActive());
339 1 : }
340 :
341 1 : TEST(NuclearDecay, superheavy_beyondZmaxSkippedSilently) {
342 : // Nuclei with Z > NUCLEAR_ZMAX (82) must be caught by the hard bounds guard
343 : // in process() before any array access, leaving the candidate unchanged.
344 1 : NuclearDecay d;
345 1 : Candidate c(nucleusId(209, 83), 80 * EeV); // Bi-209, Z=83
346 1 : c.setCurrentStep(100 * Mpc);
347 1 : int id_before = c.current.getId();
348 1 : double E_before = c.current.getEnergy();
349 1 : d.process(&c);
350 1 : EXPECT_EQ(id_before, c.current.getId());
351 1 : EXPECT_DOUBLE_EQ(E_before, c.current.getEnergy());
352 1 : EXPECT_TRUE(c.isActive());
353 1 : }
354 :
355 :
356 1 : TEST(NuclearDecay, superheavy_unstableHeavyNucleusDecays) {
357 : // Co-56 (Z=27, N=29) beta+ decays to stable Fe-56 (Z=26, N=30).
358 : // tau_rest ~ 13.9 Mdays; at 1 EeV the decay length is ~2 Mpc so it
359 : // decays well within a 100 Mpc step.
360 : // This test is skipped when the extended nuclear_decay.txt (Z<=82) is not
361 : // loaded; it can stochastically fail at energies where the decay length
362 : // exceeds the step, but the chosen energy makes that vanishingly unlikely.
363 1 : NuclearDecay d(true, false, false); // enable electrons so positron is produced
364 1 : int co56 = nucleusId(56, 27);
365 1 : if (d.meanFreePath(co56, 1e7) == std::numeric_limits<double>::max())
366 : return; // extended file not deployed
367 :
368 1 : Candidate c(co56, 1e18 * eV);
369 1 : c.setCurrentStep(100 * Mpc);
370 1 : double gamma = c.current.getLorentzFactor();
371 1 : d.process(&c);
372 :
373 : // Fe-56 is the only stable beta+ daughter of Co-56
374 1 : EXPECT_EQ(nucleusId(56, 26), c.current.getId());
375 : // Lorentz factor conserved across beta+ decay (nuclear recoil negligible)
376 1 : EXPECT_DOUBLE_EQ(gamma, c.current.getLorentzFactor());
377 : // Positron secondary from beta+ decay
378 1 : EXPECT_GE(c.secondaries.size(), 1);
379 1 : }
380 :
381 : // PhotoDisintegration --------------------------------------------------------
382 1 : TEST(PhotoDisintegration, allBackgrounds) {
383 : // Test if interaction data files are loaded.
384 1 : ref_ptr<PhotonField> cmb = new CMB();
385 1 : PhotoDisintegration pd(cmb);
386 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
387 1 : pd.setPhotonField(irb);
388 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
389 1 : pd.setPhotonField(urb);
390 1 : irb = new IRB_Stecker05();
391 1 : pd.setPhotonField(irb);
392 1 : irb = new IRB_Franceschini08();
393 1 : pd.setPhotonField(irb);
394 1 : irb = new IRB_Finke10();
395 1 : pd.setPhotonField(irb);
396 1 : irb = new IRB_Dominguez11();
397 1 : pd.setPhotonField(irb);
398 1 : irb = new IRB_Gilmore12();
399 1 : pd.setPhotonField(irb);
400 1 : irb = new IRB_Stecker16_upper();
401 1 : pd.setPhotonField(irb);
402 1 : irb = new IRB_Stecker16_lower();
403 1 : pd.setPhotonField(irb);
404 1 : irb = new IRB_Finke22();
405 1 : pd.setPhotonField(irb);
406 1 : urb = new URB_Nitu21();
407 2 : pd.setPhotonField(urb);
408 2 : }
409 :
410 1 : TEST(PhotoDisintegration, carbon) {
411 : // Test if a 100 EeV C-12 nucleus photo-disintegrates (at least once) over a distance of 1 Gpc.
412 : // This test can stochastically fail.
413 1 : ref_ptr<PhotonField> cmb = new CMB();
414 1 : PhotoDisintegration pd(cmb);
415 1 : Candidate c;
416 1 : int id = nucleusId(12, 6);
417 1 : c.current.setId(id);
418 1 : c.current.setEnergy(100 * EeV);
419 1 : c.setCurrentStep(1000 * Mpc);
420 1 : pd.process(&c);
421 :
422 1 : EXPECT_TRUE(c.current.getEnergy() < 100 * EeV);
423 : // energy loss
424 1 : EXPECT_TRUE(c.secondaries.size() > 0);
425 : // secondaries produced
426 :
427 1 : double E = c.current.getEnergy();
428 1 : id = c.current.getId();
429 1 : int A = massNumber(id);
430 1 : int Z = chargeNumber(id);
431 :
432 4 : for (int i = 0; i < c.secondaries.size(); i++) {
433 3 : E += (*c.secondaries[i]).current.getEnergy();
434 3 : id = (*c.secondaries[i]).current.getId();
435 3 : A += massNumber(id);
436 3 : Z += chargeNumber(id);
437 : }
438 1 : EXPECT_EQ(12, A);
439 : // nucleon number conserved
440 1 : EXPECT_EQ(6, Z);
441 : // proton number conserved
442 1 : EXPECT_DOUBLE_EQ(100 * EeV, E);
443 : // energy conserved
444 2 : }
445 :
446 1 : TEST(PhotoDisintegration, iron) {
447 : // Test if a 200 EeV Fe-56 nucleus photo-disintegrates (at least once) over a distance of 1 Gpc.
448 : // This test can stochastically fail.
449 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
450 1 : PhotoDisintegration pd(irb);
451 1 : Candidate c;
452 1 : int id = nucleusId(56, 26);
453 1 : c.current.setId(id);
454 1 : c.current.setEnergy(200 * EeV);
455 1 : c.setCurrentStep(1000 * Mpc);
456 1 : pd.process(&c);
457 :
458 : // expect energy loss
459 1 : EXPECT_TRUE(c.current.getEnergy() < 200 * EeV);
460 :
461 : // expect secondaries produced
462 1 : EXPECT_TRUE(c.secondaries.size() > 0);
463 :
464 1 : double E = c.current.getEnergy();
465 1 : id = c.current.getId();
466 1 : int A = massNumber(id);
467 1 : int Z = chargeNumber(id);
468 :
469 43 : for (int i = 0; i < c.secondaries.size(); i++) {
470 42 : E += (*c.secondaries[i]).current.getEnergy();
471 42 : id = (*c.secondaries[i]).current.getId();
472 42 : A += massNumber(id);
473 42 : Z += chargeNumber(id);
474 : }
475 :
476 : // nucleon number conserved
477 1 : EXPECT_EQ(56, A);
478 :
479 : // proton number conserved (no decay active)
480 1 : EXPECT_EQ(26, Z);
481 :
482 : // energy conserved
483 1 : EXPECT_DOUBLE_EQ(200 * EeV, E);
484 2 : }
485 :
486 1 : TEST(PhotoDisintegration, thisIsNotNucleonic) {
487 : // Test that nothing happens to an electron.
488 1 : ref_ptr<PhotonField> cmb = new CMB();
489 1 : PhotoDisintegration pd(cmb);
490 1 : Candidate c;
491 1 : c.setCurrentStep(1 * Mpc);
492 1 : c.current.setId(11); // electron
493 1 : c.current.setEnergy(10 * EeV);
494 1 : pd.process(&c);
495 1 : EXPECT_EQ(11, c.current.getId());
496 1 : EXPECT_EQ(10 * EeV, c.current.getEnergy());
497 2 : }
498 :
499 1 : TEST(PhotoDisintegration, limitNextStep) {
500 : // Test if the interaction limits the next propagation step.
501 1 : ref_ptr<PhotonField> cmb = new CMB();
502 1 : PhotoDisintegration pd(cmb);
503 1 : Candidate c;
504 1 : c.setNextStep(std::numeric_limits<double>::max());
505 1 : c.current.setId(nucleusId(4, 2));
506 1 : c.current.setEnergy(200 * EeV);
507 1 : pd.process(&c);
508 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
509 2 : }
510 :
511 1 : TEST(PhotoDisintegration, allIsotopes) {
512 : // Test if all isotopes are handled.
513 1 : ref_ptr<PhotonField> cmb = new CMB();
514 1 : PhotoDisintegration pd1(cmb);
515 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
516 1 : PhotoDisintegration pd2(irb);
517 1 : Candidate c;
518 1 : c.setCurrentStep(10 * Mpc);
519 :
520 27 : for (int Z = 1; Z <= 26; Z++) {
521 806 : for (int N = 1; N <= 30; N++) {
522 :
523 780 : c.current.setId(nucleusId(Z + N, Z));
524 780 : c.current.setEnergy(80 * EeV);
525 780 : pd1.process(&c);
526 :
527 780 : c.current.setId(nucleusId(Z + N, Z));
528 780 : c.current.setEnergy(80 * EeV);
529 780 : pd2.process(&c);
530 : }
531 : }
532 3 : }
533 :
534 1 : TEST(Photodisintegration, updateParticleParentProperties) { // Issue: #204
535 1 : ref_ptr<PhotonField> cmb = new CMB();
536 1 : PhotoDisintegration pd(cmb);
537 :
538 1 : Candidate c(nucleusId(56,26), 500 * EeV, Vector3d(1 * Mpc, 0, 0));
539 :
540 1 : pd.performInteraction(&c, 1);
541 : // the candidates parent is the original particle
542 1 : EXPECT_EQ(c.created.getId(), nucleusId(56,26));
543 :
544 1 : pd.performInteraction(&c, 1);
545 : // now it has to be changed
546 1 : EXPECT_NE(c.created.getId(), nucleusId(56,26));
547 2 : }
548 :
549 1 : TEST(PhotoDisintegration, interactionTag) {
550 1 : PhotoDisintegration pd(new CMB());
551 :
552 : // test default interactionTag
553 1 : EXPECT_TRUE(pd.getInteractionTag() == "PD");
554 :
555 : // test secondary tag
556 1 : pd.setHavePhotons(true);
557 2 : Candidate c(nucleusId(56,26), 500 * EeV);
558 1 : c.setCurrentStep(1 * Gpc);
559 1 : pd.process(&c);
560 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "PD");
561 :
562 : // test custom tag
563 1 : pd.setInteractionTag("myTag");
564 1 : EXPECT_TRUE(pd.getInteractionTag() == "myTag");
565 1 : }
566 :
567 : // PhotoDisintegration - Superheavy extension ---------------------------------
568 :
569 1 : TEST(PhotoDisintegration, superheavy_standardTablesSkipHeavyNuclei) {
570 : // With standard tables (superheavy=false), nuclei with Z > 26 have no
571 : // entries in the TALYS 1.8 rate files. The empty-rate guard in process()
572 : // must silently skip them, leaving the candidate unchanged.
573 1 : ref_ptr<PhotonField> cmb = new CMB();
574 1 : PhotoDisintegration pd(cmb, false, 0.1, false);
575 1 : Candidate c;
576 1 : c.current.setId(nucleusId(58, 28)); // Ni-58, Z=28
577 1 : c.current.setEnergy(200 * EeV);
578 1 : c.setCurrentStep(1000 * Mpc);
579 1 : int id_before = c.current.getId();
580 1 : double E_before = c.current.getEnergy();
581 1 : pd.process(&c);
582 1 : EXPECT_EQ(id_before, c.current.getId());
583 1 : EXPECT_DOUBLE_EQ(E_before, c.current.getEnergy());
584 2 : }
585 :
586 1 : TEST(PhotoDisintegration, superheavy_boundsGuardZmax) {
587 : // Nuclei with Z > NUCLEAR_ZMAX (82) must be caught by the hard bounds
588 : // guard in both process() and lossLength() before any array access.
589 : // The candidate must be left completely unchanged by process().
590 1 : ref_ptr<PhotonField> cmb = new CMB();
591 1 : PhotoDisintegration pd(cmb);
592 : // Bi-209: Z=83 > NUCLEAR_ZMAX=82
593 1 : Candidate c;
594 1 : c.current.setId(nucleusId(209, 83));
595 1 : c.current.setEnergy(200 * EeV);
596 1 : c.setCurrentStep(1000 * Mpc);
597 1 : int id_before = c.current.getId();
598 1 : double E_before = c.current.getEnergy();
599 1 : pd.process(&c);
600 1 : EXPECT_EQ(id_before, c.current.getId());
601 1 : EXPECT_DOUBLE_EQ(E_before, c.current.getEnergy());
602 2 : }
603 :
604 1 : TEST(PhotoDisintegration, superheavy_boundsGuardNmax) {
605 : // Nuclei with N > NUCLEAR_NMAX (132) must be caught by the hard bounds
606 : // guard in process() before any array access.
607 1 : ref_ptr<PhotonField> cmb = new CMB();
608 1 : PhotoDisintegration pd(cmb);
609 : // H-134: Z=1, N=133 > NUCLEAR_NMAX=132
610 1 : Candidate c;
611 1 : c.current.setId(nucleusId(134, 1));
612 1 : c.current.setEnergy(200 * EeV);
613 1 : c.setCurrentStep(1000 * Mpc);
614 1 : int id_before = c.current.getId();
615 1 : double E_before = c.current.getEnergy();
616 1 : pd.process(&c);
617 1 : EXPECT_EQ(id_before, c.current.getId());
618 1 : EXPECT_DOUBLE_EQ(E_before, c.current.getEnergy());
619 2 : }
620 :
621 :
622 1 : TEST(PhotoDisintegration, superheavy_lead) {
623 : // Test if a Pb-208 nucleus photo-disintegrates at least once over 1 Gpc
624 : // using the superheavy CMB tables. A, Z, and energy conservation are
625 : // verified across all interactions in the step.
626 : // This test can stochastically fail.
627 1 : ref_ptr<PhotonField> cmb = new CMB();
628 1 : PhotoDisintegration pd(cmb, false, 0.1, true);
629 1 : Candidate c;
630 1 : int id = nucleusId(208, 82);
631 1 : c.current.setId(id);
632 : // Energy chosen so log10(gamma) ~ 10 for Pb-208, above the CMB threshold.
633 1 : c.current.setEnergy(2000 * EeV);
634 1 : c.setCurrentStep(1000 * Mpc);
635 1 : pd.process(&c);
636 :
637 1 : EXPECT_LT(c.current.getEnergy(), 2000 * EeV);
638 1 : EXPECT_GT(c.secondaries.size(), 0);
639 :
640 1 : double E = c.current.getEnergy();
641 1 : id = c.current.getId();
642 1 : int A = massNumber(id);
643 1 : int Z = chargeNumber(id);
644 14 : for (size_t i = 0; i < c.secondaries.size(); i++) {
645 13 : E += (*c.secondaries[i]).current.getEnergy();
646 13 : id = (*c.secondaries[i]).current.getId();
647 13 : A += massNumber(id);
648 13 : Z += chargeNumber(id);
649 : }
650 : // nucleon number conserved
651 1 : EXPECT_EQ(208, A);
652 : // proton number conserved
653 1 : EXPECT_EQ(82, Z);
654 : // energy conserved (EXPECT_NEAR: multiple interactions accumulate ~1e-12 J rounding)
655 1 : EXPECT_NEAR(2000 * EeV, E, 1e-9);
656 2 : }
657 :
658 1 : TEST(PhotoDisintegration, superheavy_allHeavyIsotopes) {
659 : // With superheavy=true, processing any nucleus in Z=27..NUCLEAR_ZMAX must
660 : // not crash, even for isotopes whose slots are empty in the rate table.
661 1 : ref_ptr<PhotonField> cmb = new CMB();
662 1 : PhotoDisintegration pd(cmb, false, 0.1, true);
663 1 : Candidate c;
664 1 : c.setCurrentStep(10 * Mpc);
665 57 : for (int Z = 27; Z <= NUCLEAR_ZMAX; Z++) {
666 1736 : for (int N = 1; N <= 30; N++) {
667 1680 : c.current.setId(nucleusId(Z + N, Z));
668 1680 : c.current.setEnergy(80 * EeV);
669 1680 : pd.process(&c);
670 : }
671 : }
672 2 : }
673 :
674 1 : TEST(PhotoDisintegration, superheavy_setPhotonFieldReloads) {
675 : // setPhotonField(field, true/false) must swap between superheavy and standard
676 : // tables. Pb-208 (only in superheavy tables) is used as the discriminating
677 : // nucleus: with standard tables it has no rate so it must not limit the next
678 : // step; with superheavy tables it has a rate so it must limit the next step.
679 1 : ref_ptr<PhotonField> cmb = new CMB();
680 1 : PhotoDisintegration pd(cmb); // start with standard tables
681 :
682 : // Standard: Pb-208 has no rate — step must not be limited
683 1 : Candidate c_std(nucleusId(208, 82), 2000 * EeV);
684 1 : c_std.setNextStep(std::numeric_limits<double>::max());
685 1 : pd.process(&c_std);
686 1 : EXPECT_DOUBLE_EQ(c_std.getNextStep(), std::numeric_limits<double>::max());
687 :
688 : // Switch to superheavy: Pb-208 now has a rate — step must be limited
689 1 : pd.setPhotonField(cmb, true);
690 2 : Candidate c_shv(nucleusId(208, 82), 2000 * EeV);
691 1 : c_shv.setNextStep(std::numeric_limits<double>::max());
692 1 : pd.process(&c_shv);
693 1 : EXPECT_LT(c_shv.getNextStep(), std::numeric_limits<double>::max());
694 :
695 : // Switch back to standard: step no longer limited for Pb-208
696 1 : pd.setPhotonField(cmb, false);
697 2 : Candidate c_back(nucleusId(208, 82), 2000 * EeV);
698 1 : c_back.setNextStep(std::numeric_limits<double>::max());
699 1 : pd.process(&c_back);
700 1 : EXPECT_DOUBLE_EQ(c_back.getNextStep(), std::numeric_limits<double>::max());
701 2 : }
702 :
703 1 : TEST(PhotoDisintegration, superheavy_overlapRatesConsistent) {
704 : // Nuclei common to both the standard rate tables and the superheavy
705 : // tables must yield consistent energy loss lengths (within 1%).
706 : // The cross sections branches for A<=12 are using the same existing
707 : // data since CRPropa 2.0
708 1 : ref_ptr<PhotonField> cmb = new CMB();
709 1 : ref_ptr<PhotonField> irb = new IRB_Gilmore12();
710 2 : PhotoDisintegration pd_cmb_std(cmb);
711 2 : PhotoDisintegration pd_cmb_shv(cmb, false, 0.1, true);
712 2 : PhotoDisintegration pd_irb_std(irb);
713 1 : PhotoDisintegration pd_irb_shv(irb, false, 0.1, true);
714 :
715 : // Excluded from this comparison:
716 : // A<12 — cross section tables unchanged;
717 : // N-14 — intrinsic ~1.77% cross-section difference reflected in
718 : // the rate tables itself; outlier as all other nuclei agree
719 : // to within <0.025%.
720 1 : struct Nucleus { int A, Z; } nuclei[] = {
721 : {12, 6},
722 : {16, 8},
723 : {20, 10},
724 : {24, 12},
725 : {28, 14},
726 : {40, 20},
727 : {56, 26},
728 : };
729 :
730 : // Sweep log10(gamma) = 6..11 using the sum of CMB + IRB rates.
731 : // Above lg=11 both datasets begin to diverge due to extrapolation, so the
732 : // range is capped there, suitable for most cosmic ray applications.
733 : // Tolerance is 0.1%: actual deviations for the listed nuclei are <0.025%,
734 : // giving a 4x margin.
735 : const int nSteps = 126;
736 : const double lgMin = 6.0, lgMax = 11.0;
737 : const double llMax = std::numeric_limits<double>::max();
738 :
739 8 : for (auto &nuc : nuclei) {
740 7 : int id = nucleusId(nuc.A, nuc.Z);
741 : double sumRatio = 0;
742 : int nValid = 0;
743 889 : for (int i = 0; i < nSteps; i++) {
744 882 : double gamma = pow(10, lgMin + (lgMax - lgMin) * i / (nSteps - 1));
745 882 : double ll_cmb_std = pd_cmb_std.lossLength(id, gamma);
746 882 : double ll_cmb_shv = pd_cmb_shv.lossLength(id, gamma);
747 882 : double ll_irb_std = pd_irb_std.lossLength(id, gamma);
748 882 : double ll_irb_shv = pd_irb_shv.lossLength(id, gamma);
749 882 : double rate_std = (ll_cmb_std < llMax ? 1/ll_cmb_std : 0)
750 882 : + (ll_irb_std < llMax ? 1/ll_irb_std : 0);
751 882 : double rate_shv = (ll_cmb_shv < llMax ? 1/ll_cmb_shv : 0)
752 882 : + (ll_irb_shv < llMax ? 1/ll_irb_shv : 0);
753 882 : if (rate_std == 0 || rate_shv == 0) continue;
754 875 : sumRatio += rate_std / rate_shv;
755 875 : nValid++;
756 : }
757 7 : if (nValid == 0) continue;
758 7 : double meanRatio = sumRatio / nValid;
759 7 : EXPECT_GT(meanRatio, 0.999) << "A=" << nuc.A << " Z=" << nuc.Z
760 : << " mean ratio=" << meanRatio << " over " << nValid << " boosts";
761 7 : EXPECT_LT(meanRatio, 1.001) << "A=" << nuc.A << " Z=" << nuc.Z
762 : << " mean ratio=" << meanRatio << " over " << nValid << " boosts";
763 : }
764 2 : }
765 :
766 : // ElasticScattering ----------------------------------------------------------
767 1 : TEST(ElasticScattering, allBackgrounds) {
768 : // Test if interaction data files are loaded.
769 1 : ref_ptr<PhotonField> cmb = new CMB();
770 1 : ElasticScattering scattering(cmb);
771 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
772 1 : scattering.setPhotonField(irb);
773 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
774 2 : scattering.setPhotonField(urb);
775 2 : }
776 :
777 1 : TEST(ElasticScattering, secondaries) {
778 : // Test the creation of cosmic ray photons.
779 : // This test can stochastically fail.
780 1 : ref_ptr<PhotonField> cmb = new CMB();
781 1 : ElasticScattering scattering(cmb);
782 1 : Candidate c;
783 1 : int id = nucleusId(12, 6);
784 1 : c.current.setId(id);
785 1 : c.current.setEnergy(200 * EeV);
786 1 : c.setCurrentStep(400 * Mpc);
787 1 : scattering.process(&c);
788 :
789 1 : EXPECT_GT(c.secondaries.size(), 0);
790 :
791 10 : for (int i = 0; i < c.secondaries.size(); i++) {
792 9 : int id = (*c.secondaries[i]).current.getId();
793 9 : EXPECT_EQ(id, 22);
794 9 : double energy = (*c.secondaries[i]).current.getEnergy();
795 9 : EXPECT_GT(energy, 0);
796 9 : EXPECT_LT(energy, 200 * EeV);
797 : }
798 2 : }
799 :
800 : // PhotoPionProduction --------------------------------------------------------
801 1 : TEST(PhotoPionProduction, allBackgrounds) {
802 : // Test if all interaction data files can be loaded.
803 1 : ref_ptr<PhotonField> cmb = new CMB();
804 1 : PhotoPionProduction ppp(cmb);
805 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
806 1 : ppp.setPhotonField(irb);
807 1 : irb = new IRB_Stecker05();
808 1 : ppp.setPhotonField(irb);
809 1 : irb = new IRB_Franceschini08();
810 1 : ppp.setPhotonField(irb);
811 1 : irb = new IRB_Finke10();
812 1 : ppp.setPhotonField(irb);
813 1 : irb = new IRB_Dominguez11();
814 1 : ppp.setPhotonField(irb);
815 1 : irb = new IRB_Gilmore12();
816 1 : ppp.setPhotonField(irb);
817 1 : irb = new IRB_Stecker16_upper();
818 1 : ppp.setPhotonField(irb);
819 1 : irb = new IRB_Stecker16_lower();
820 1 : ppp.setPhotonField(irb);
821 1 : irb = new IRB_Finke22();
822 1 : ppp.setPhotonField(irb);
823 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
824 1 : ppp.setPhotonField(urb);
825 1 : urb = new URB_Nitu21();
826 2 : ppp.setPhotonField(urb);
827 2 : }
828 :
829 1 : TEST(PhotoPionProduction, proton) {
830 : // Test photopion interaction for 100 EeV proton.
831 : // This test can stochastically fail.
832 1 : ref_ptr<PhotonField> cmb = new CMB();
833 1 : PhotoPionProduction ppp(cmb);
834 1 : Candidate c(nucleusId(1, 1), 100 * EeV);
835 1 : c.setCurrentStep(1000 * Mpc);
836 1 : ppp.process(&c);
837 :
838 : // expect energy loss
839 1 : EXPECT_LT(c.current.getEnergy(), 100. * EeV);
840 :
841 : // expect nucleon number conservation
842 1 : EXPECT_EQ(1, massNumber(c.current.getId()));
843 :
844 : // expect no (nucleonic) secondaries
845 1 : EXPECT_EQ(0, c.secondaries.size());
846 2 : }
847 :
848 1 : TEST(PhotoPionProduction, helium) {
849 : // Test photo-pion interaction for 400 EeV He nucleus.
850 : // This test can stochastically fail.
851 1 : ref_ptr<PhotonField> cmb = new CMB();
852 1 : PhotoPionProduction ppp(cmb);
853 1 : Candidate c;
854 1 : c.current.setId(nucleusId(4, 2));
855 1 : c.current.setEnergy(400. * EeV);
856 1 : c.setCurrentStep(1000 * Mpc);
857 1 : ppp.process(&c);
858 1 : EXPECT_LT(c.current.getEnergy(), 400. * EeV);
859 1 : int id = c.current.getId();
860 1 : EXPECT_TRUE(massNumber(id) < 4);
861 1 : EXPECT_TRUE(c.secondaries.size() > 0);
862 2 : }
863 :
864 1 : TEST(PhotoPionProduction, thisIsNotNucleonic) {
865 : // Test if nothing happens to an electron.
866 1 : ref_ptr<PhotonField> cmb = new CMB();
867 1 : PhotoPionProduction ppp(cmb);
868 1 : Candidate c;
869 1 : c.current.setId(11); // electron
870 1 : c.current.setEnergy(10 * EeV);
871 1 : c.setCurrentStep(100 * Mpc);
872 1 : ppp.process(&c);
873 1 : EXPECT_EQ(11, c.current.getId());
874 1 : EXPECT_EQ(10 * EeV, c.current.getEnergy());
875 2 : }
876 :
877 1 : TEST(PhotoPionProduction, limitNextStep) {
878 : // Test if the interaction limits the next propagation step.
879 1 : ref_ptr<PhotonField> cmb = new CMB();
880 1 : PhotoPionProduction ppp(cmb);
881 1 : Candidate c(nucleusId(1, 1), 200 * EeV);
882 1 : c.setNextStep(std::numeric_limits<double>::max());
883 1 : ppp.process(&c);
884 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
885 2 : }
886 :
887 1 : TEST(PhotoPionProduction, secondaries) {
888 : // Test photo-pion interaction for 100 EeV proton.
889 : // This test can stochastically fail.
890 1 : ref_ptr<PhotonField> cmb = new CMB();
891 1 : PhotoPionProduction ppp(cmb, true, true, true);
892 1 : Candidate c(nucleusId(1, 1), 100 * EeV);
893 1 : c.setCurrentStep(1000 * Mpc);
894 1 : ppp.process(&c);
895 : // there should be secondaries
896 1 : EXPECT_GT(c.secondaries.size(), 1);
897 2 : }
898 :
899 1 : TEST(PhotoPionProduction, sampling) {
900 : // Specific test of photon sampling of photo-pion production
901 : // by testing the calculated pEpsMax for CMB(), also indirectly
902 : // testing epsMinInteraction and logSampling (default).
903 1 : ref_ptr<PhotonField> cmb = new CMB(); //create CMB instance
904 : double energy = 1.e10; //1e10 GeV
905 : bool onProton = true; //proton
906 : double z = 0; //no redshift
907 1 : PhotoPionProduction ppp(cmb, true, true, true);
908 1 : double correctionFactor = ppp.getCorrectionFactor(); //get current correctionFactor
909 2 : double epsMin = std::max(cmb -> getMinimumPhotonEnergy(z) / eV, 0.00710614); // 0.00710614 = epsMinInteraction(onProton,energy)
910 1 : double epsMax = cmb -> getMaximumPhotonEnergy(z) / eV;
911 1 : double pEpsMax = ppp.probEpsMax(onProton, energy, z, epsMin, epsMax) / correctionFactor;
912 1 : EXPECT_DOUBLE_EQ(pEpsMax,132673934934.922);
913 2 : }
914 :
915 1 : TEST(PhotoPionProduction, interactionTag) {
916 1 : PhotoPionProduction ppp(new CMB());
917 :
918 : // test default interactionTag
919 1 : EXPECT_TRUE(ppp.getInteractionTag() == "PPP");
920 :
921 : // test secondary tag
922 1 : ppp.setHavePhotons(true);
923 2 : Candidate c(nucleusId(1,1), 100 * EeV);
924 11 : for(int i = 0; i <10; i++)
925 10 : ppp.performInteraction(&c, true);
926 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "PPP");
927 :
928 : // test custom interactionTag
929 1 : ppp.setInteractionTag("myTag");
930 1 : EXPECT_TRUE(ppp.getInteractionTag() == "myTag");
931 1 : }
932 :
933 : // Redshift -------------------------------------------------------------------
934 2 : TEST(Redshift, simpleTest) {
935 : // Test if redshift is decreased and adiabatic energy loss is applied.
936 : Redshift redshift;
937 :
938 1 : Candidate c;
939 1 : c.setRedshift(0.024);
940 1 : c.current.setEnergy(100 * EeV);
941 1 : c.setCurrentStep(1 * Mpc);
942 :
943 1 : redshift.process(&c);
944 1 : EXPECT_GT(0.024, c.getRedshift()); // expect redshift decrease
945 1 : EXPECT_GT(100, c.current.getEnergy() / EeV); // expect energy loss
946 2 : }
947 :
948 2 : TEST(Redshift, limitRedshiftDecrease) {
949 : // Test if the redshift decrease is limited to z_updated >= 0.
950 : Redshift redshift;
951 :
952 1 : Candidate c;
953 1 : c.setRedshift(0.024); // roughly corresponds to 100 Mpc
954 1 : c.setCurrentStep(150 * Mpc);
955 :
956 1 : redshift.process(&c);
957 1 : EXPECT_DOUBLE_EQ(0, c.getRedshift());
958 2 : }
959 :
960 : // EMPairProduction -----------------------------------------------------------
961 1 : TEST(EMPairProduction, allBackgrounds) {
962 : // Test if interaction data files are loaded.
963 1 : ref_ptr<PhotonField> cmb = new CMB();
964 2 : EMPairProduction em(cmb);
965 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
966 1 : em.setPhotonField(ebl);
967 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
968 1 : em.setPhotonField(urb);
969 :
970 1 : ebl = new IRB_Stecker05();
971 1 : em.setPhotonField(ebl);
972 1 : ebl = new IRB_Franceschini08();
973 1 : em.setPhotonField(ebl);
974 1 : ebl = new IRB_Finke10();
975 1 : em.setPhotonField(ebl);
976 1 : ebl = new IRB_Dominguez11();
977 1 : em.setPhotonField(ebl);
978 1 : ebl = new IRB_Gilmore12();
979 1 : em.setPhotonField(ebl);
980 1 : ebl = new IRB_Stecker16_upper();
981 1 : em.setPhotonField(ebl);
982 1 : ebl = new IRB_Stecker16_lower();
983 1 : em.setPhotonField(ebl);
984 1 : ebl = new IRB_Finke22();
985 1 : em.setPhotonField(ebl);
986 1 : urb = new URB_Fixsen11();
987 1 : em.setPhotonField(urb);
988 1 : urb = new URB_Nitu21();
989 2 : em.setPhotonField(urb);
990 2 : }
991 :
992 1 : TEST(EMPairProduction, limitNextStep) {
993 : // Test if the interaction limits the next propagation step.
994 1 : ref_ptr<PhotonField> cmb = new CMB();
995 2 : EMPairProduction m(cmb);
996 1 : Candidate c(22, 1E17 * eV);
997 1 : c.setNextStep(std::numeric_limits<double>::max());
998 1 : m.process(&c);
999 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
1000 2 : }
1001 :
1002 1 : TEST(EMPairProduction, secondaries) {
1003 : // Test if secondaries are correctly produced.
1004 1 : ref_ptr<PhotonField> cmb = new CMB();
1005 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
1006 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
1007 2 : EMPairProduction m(cmb);
1008 1 : m.setHaveElectrons(true);
1009 1 : m.setThinning(0.);
1010 :
1011 : std::vector<ref_ptr<PhotonField>> fields;
1012 1 : fields.push_back(cmb);
1013 1 : fields.push_back(irb);
1014 1 : fields.push_back(urb);
1015 :
1016 : // loop over photon backgrounds
1017 4 : for (int f = 0; f < fields.size(); f++) {
1018 3 : m.setPhotonField(fields[f]);
1019 423 : for (int i = 0; i < 140; i++) { // loop over energies Ep = (1e10 - 1e23) eV
1020 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
1021 420 : Candidate c(22, Ep);
1022 420 : c.setCurrentStep(1e4 * Mpc);
1023 420 : m.process(&c);
1024 :
1025 : // pass if no interaction has ocurred (no tabulated rates)
1026 420 : if (c.isActive())
1027 : continue;
1028 :
1029 : // expect 2 secondaries
1030 245 : EXPECT_EQ(c.secondaries.size(), 2);
1031 :
1032 : // expect electron / positron with energies 0 < E < Ephoton
1033 : double Etot = 0;
1034 735 : for (int j = 0; j < c.secondaries.size(); j++) {
1035 490 : Candidate s = *c.secondaries[j];
1036 490 : EXPECT_EQ(abs(s.current.getId()), 11);
1037 490 : EXPECT_GT(s.current.getEnergy(), 0);
1038 490 : EXPECT_LT(s.current.getEnergy(), Ep);
1039 490 : Etot += s.current.getEnergy();
1040 490 : }
1041 :
1042 : // test energy conservation
1043 245 : EXPECT_DOUBLE_EQ(Ep, Etot);
1044 420 : }
1045 : }
1046 2 : }
1047 :
1048 2 : TEST(EMPairProduction, interactionTag) {
1049 2 : EMPairProduction m(new CMB());
1050 :
1051 : // test default interactionTag
1052 1 : EXPECT_TRUE(m.getInteractionTag() == "EMPP");
1053 :
1054 : // test secondary tag
1055 1 : m.setHaveElectrons(true);
1056 2 : Candidate c(22, 1 * EeV);
1057 1 : m.performInteraction(&c);
1058 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMPP");
1059 :
1060 : // test custom tag
1061 1 : m.setInteractionTag("myTag");
1062 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
1063 1 : }
1064 :
1065 : // EMDoublePairProduction -----------------------------------------------------
1066 1 : TEST(EMDoublePairProduction, allBackgrounds) {
1067 : // Test if interaction data files are loaded.
1068 1 : ref_ptr<PhotonField> cmb = new CMB();
1069 2 : EMDoublePairProduction em(cmb);
1070 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
1071 1 : em.setPhotonField(ebl);
1072 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
1073 1 : em.setPhotonField(urb);
1074 :
1075 1 : ebl = new IRB_Stecker05();
1076 1 : em.setPhotonField(ebl);
1077 1 : ebl = new IRB_Franceschini08();
1078 1 : em.setPhotonField(ebl);
1079 1 : ebl = new IRB_Finke10();
1080 1 : em.setPhotonField(ebl);
1081 1 : ebl = new IRB_Dominguez11();
1082 1 : em.setPhotonField(ebl);
1083 1 : ebl = new IRB_Gilmore12();
1084 1 : em.setPhotonField(ebl);
1085 1 : ebl = new IRB_Stecker16_upper();
1086 1 : em.setPhotonField(ebl);
1087 1 : ebl = new IRB_Stecker16_lower();
1088 1 : em.setPhotonField(ebl);
1089 1 : ebl = new IRB_Finke22();
1090 1 : em.setPhotonField(ebl);
1091 1 : urb = new URB_Fixsen11();
1092 1 : em.setPhotonField(urb);
1093 1 : urb = new URB_Nitu21();
1094 2 : em.setPhotonField(urb);
1095 2 : }
1096 :
1097 1 : TEST(EMDoublePairProduction, limitNextStep) {
1098 : // Test if the interaction limits the next propagation step.
1099 1 : ref_ptr<PhotonField> cmb = new CMB();
1100 2 : EMDoublePairProduction m(cmb);
1101 1 : Candidate c(22, 1E17 * eV);
1102 1 : c.setNextStep(std::numeric_limits<double>::max());
1103 1 : m.process(&c);
1104 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
1105 2 : }
1106 :
1107 1 : TEST(EMDoublePairProduction, secondaries) {
1108 : // Test if secondaries are correctly produced.
1109 1 : ref_ptr<PhotonField> cmb = new CMB();
1110 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
1111 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
1112 2 : EMPairProduction m(cmb);
1113 1 : m.setHaveElectrons(true);
1114 1 : m.setThinning(0.);
1115 :
1116 : std::vector<ref_ptr<PhotonField>> fields;
1117 1 : fields.push_back(cmb);
1118 1 : fields.push_back(irb);
1119 1 : fields.push_back(urb);
1120 :
1121 : // loop over photon backgrounds
1122 4 : for (int f = 0; f < fields.size(); f++) {
1123 3 : m.setPhotonField(fields[f]);
1124 :
1125 : // loop over energies Ep = (1e9 - 1e23) eV
1126 423 : for (int i = 0; i < 140; i++) {
1127 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
1128 420 : Candidate c(22, Ep);
1129 420 : c.setCurrentStep(1e4 * Mpc); // use lower value so that the test can run faster
1130 420 : m.process(&c);
1131 :
1132 : // pass if no interaction has occured (no tabulated rates)
1133 420 : if (c.isActive())
1134 : continue;
1135 :
1136 : // expect 2 secondaries (only one pair is considered)
1137 252 : EXPECT_EQ(c.secondaries.size(), 2);
1138 :
1139 : // expect electron / positron with energies 0 < E < Ephoton
1140 : double Etot = 0;
1141 756 : for (int j = 0; j < c.secondaries.size(); j++) {
1142 504 : Candidate s = *c.secondaries[j];
1143 504 : EXPECT_EQ(abs(s.current.getId()), 11);
1144 504 : EXPECT_GT(s.current.getEnergy(), 0);
1145 504 : EXPECT_LT(s.current.getEnergy(), Ep);
1146 504 : Etot += s.current.getEnergy();
1147 504 : }
1148 :
1149 : // test energy conservation
1150 252 : EXPECT_NEAR(Ep, Etot, 1E-9);
1151 420 : }
1152 : }
1153 2 : }
1154 :
1155 2 : TEST(EMDoublePairProduction, interactionTag) {
1156 2 : EMDoublePairProduction m(new CMB());
1157 :
1158 : // test default interactionTag
1159 1 : EXPECT_TRUE(m.getInteractionTag() == "EMDP");
1160 :
1161 : // test secondary tag
1162 1 : m.setHaveElectrons(true);
1163 2 : Candidate c(22, 1 * EeV);
1164 1 : m.performInteraction(&c);
1165 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMDP");
1166 :
1167 : // test custom tag
1168 1 : m.setInteractionTag("myTag");
1169 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
1170 1 : }
1171 :
1172 : // EMTripletPairProduction ----------------------------------------------------
1173 1 : TEST(EMTripletPairProduction, allBackgrounds) {
1174 : // Test if interaction data files are loaded.
1175 1 : ref_ptr<PhotonField> cmb = new CMB();
1176 2 : EMTripletPairProduction em(cmb);
1177 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
1178 1 : em.setPhotonField(ebl);
1179 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
1180 1 : em.setPhotonField(urb);
1181 :
1182 1 : ebl = new IRB_Stecker05();
1183 1 : em.setPhotonField(ebl);
1184 1 : ebl = new IRB_Franceschini08();
1185 1 : em.setPhotonField(ebl);
1186 1 : ebl = new IRB_Finke10();
1187 1 : em.setPhotonField(ebl);
1188 1 : ebl = new IRB_Dominguez11();
1189 1 : em.setPhotonField(ebl);
1190 1 : ebl = new IRB_Gilmore12();
1191 1 : em.setPhotonField(ebl);
1192 1 : ebl = new IRB_Stecker16_upper();
1193 1 : em.setPhotonField(ebl);
1194 1 : ebl = new IRB_Stecker16_lower();
1195 1 : em.setPhotonField(ebl);
1196 1 : ebl = new IRB_Finke22();
1197 1 : em.setPhotonField(ebl);
1198 1 : urb = new URB_Fixsen11();
1199 1 : em.setPhotonField(urb);
1200 1 : urb = new URB_Nitu21();
1201 2 : em.setPhotonField(urb);
1202 2 : }
1203 :
1204 1 : TEST(EMTripletPairProduction, limitNextStep) {
1205 : // Test if the interaction limits the next propagation step.
1206 1 : ref_ptr<PhotonField> cmb = new CMB();
1207 2 : EMTripletPairProduction m(cmb);
1208 1 : Candidate c(11, 1E17 * eV);
1209 1 : c.setNextStep(std::numeric_limits<double>::max());
1210 1 : m.process(&c);
1211 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
1212 2 : }
1213 :
1214 1 : TEST(EMTripletPairProduction, secondaries) {
1215 : // Test if secondaries are correctly produced.
1216 1 : ref_ptr<PhotonField> cmb = new CMB();
1217 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
1218 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
1219 2 : EMPairProduction m(cmb);
1220 1 : m.setHaveElectrons(true);
1221 1 : m.setThinning(0.);
1222 :
1223 : std::vector<ref_ptr<PhotonField>> fields;
1224 1 : fields.push_back(cmb);
1225 1 : fields.push_back(irb);
1226 1 : fields.push_back(urb);
1227 :
1228 : // loop over photon backgrounds
1229 4 : for (int f = 0; f < fields.size(); f++) {
1230 3 : m.setPhotonField(fields[f]);
1231 :
1232 : // loop over energies Ep = (1e9 - 1e23) eV
1233 423 : for (int i = 0; i < 140; i++) {
1234 :
1235 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
1236 420 : Candidate c(11, Ep);
1237 420 : c.setCurrentStep(1e4 * Mpc); // use lower value so that the test can run faster
1238 420 : m.process(&c);
1239 :
1240 : // pass if no interaction has occured (no tabulated rates)
1241 420 : if (c.current.getEnergy() == Ep)
1242 : continue;
1243 :
1244 : // expect positive energy of primary electron
1245 0 : EXPECT_GT(c.current.getEnergy(), 0);
1246 0 : double Etot = c.current.getEnergy();
1247 :
1248 : // expect electron / positron with energies 0 < E < Ephoton
1249 0 : for (int j = 0; j < c.secondaries.size(); j++) {
1250 0 : Candidate s = *c.secondaries[j];
1251 0 : EXPECT_EQ(abs(s.current.getId()), 11);
1252 0 : EXPECT_GT(s.current.getEnergy(), 0);
1253 0 : EXPECT_LT(s.current.getEnergy(), Ep);
1254 0 : Etot += s.current.getEnergy();
1255 0 : }
1256 :
1257 : // test energy conservation
1258 0 : EXPECT_NEAR(Ep, Etot, 1e-9);
1259 420 : }
1260 : }
1261 2 : }
1262 :
1263 2 : TEST(EMTripletPairProduction, interactionTag) {
1264 2 : EMTripletPairProduction m(new CMB());
1265 :
1266 : // test default interactionTag
1267 1 : EXPECT_TRUE(m.getInteractionTag() == "EMTP");
1268 :
1269 : // test secondary tag
1270 1 : m.setHaveElectrons(true);
1271 2 : Candidate c(11, 1 * EeV);
1272 1 : m.performInteraction(&c);
1273 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMTP");
1274 :
1275 : // test custom tag
1276 1 : m.setInteractionTag("myTag");
1277 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
1278 1 : }
1279 :
1280 : // EMInverseComptonScattering -------------------------------------------------
1281 1 : TEST(EMInverseComptonScattering, allBackgrounds) {
1282 : // Test if interaction data files are loaded.
1283 1 : ref_ptr<PhotonField> cmb = new CMB();
1284 2 : EMInverseComptonScattering em(cmb);
1285 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
1286 1 : em.setPhotonField(ebl);
1287 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
1288 1 : em.setPhotonField(urb);
1289 :
1290 1 : ebl = new IRB_Stecker05();
1291 1 : em.setPhotonField(ebl);
1292 1 : ebl = new IRB_Franceschini08();
1293 1 : em.setPhotonField(ebl);
1294 1 : ebl = new IRB_Finke10();
1295 1 : em.setPhotonField(ebl);
1296 1 : ebl = new IRB_Dominguez11();
1297 1 : em.setPhotonField(ebl);
1298 1 : ebl = new IRB_Gilmore12();
1299 1 : em.setPhotonField(ebl);
1300 1 : ebl = new IRB_Stecker16_upper();
1301 1 : em.setPhotonField(ebl);
1302 1 : ebl = new IRB_Stecker16_lower();
1303 1 : em.setPhotonField(ebl);
1304 1 : ebl = new IRB_Finke22();
1305 1 : em.setPhotonField(ebl);
1306 1 : urb = new URB_Fixsen11();
1307 1 : em.setPhotonField(urb);
1308 1 : urb = new URB_Nitu21();
1309 2 : em.setPhotonField(urb);
1310 2 : }
1311 :
1312 1 : TEST(EMInverseComptonScattering, limitNextStep) {
1313 : // Test if the interaction limits the next propagation step.
1314 1 : ref_ptr<PhotonField> cmb = new CMB();
1315 2 : EMInverseComptonScattering m(cmb);
1316 1 : Candidate c(11, 1E17 * eV);
1317 1 : c.setNextStep(std::numeric_limits<double>::max());
1318 1 : m.process(&c);
1319 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
1320 2 : }
1321 :
1322 1 : TEST(EMInverseComptonScattering, secondaries) {
1323 : // Test if secondaries are correctly produced.
1324 1 : ref_ptr<PhotonField> cmb = new CMB();
1325 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
1326 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
1327 2 : EMPairProduction m(cmb);
1328 1 : m.setHaveElectrons(true);
1329 1 : m.setThinning(0.);
1330 :
1331 : std::vector<ref_ptr<PhotonField>> fields;
1332 1 : fields.push_back(cmb);
1333 1 : fields.push_back(irb);
1334 1 : fields.push_back(urb);
1335 :
1336 : // loop over photon backgrounds
1337 4 : for (int f = 0; f < fields.size(); f++) {
1338 3 : m.setPhotonField(fields[f]);
1339 :
1340 : // loop over energies Ep = (1e9 - 1e23) eV
1341 423 : for (int i = 0; i < 140; i++) {
1342 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
1343 420 : Candidate c(11, Ep);
1344 420 : c.setCurrentStep(1e3 * Mpc); // use lower value so that the test can run faster
1345 420 : m.process(&c);
1346 :
1347 : // pass if no interaction has occured (no tabulated rates)
1348 420 : if (c.current.getEnergy() == Ep)
1349 : continue;
1350 :
1351 : // expect positive energy of primary electron
1352 0 : EXPECT_GT(c.current.getEnergy(), 0);
1353 :
1354 : // expect photon with energy 0 < E < Ephoton
1355 0 : Candidate s = *c.secondaries[0];
1356 0 : EXPECT_EQ(abs(s.current.getId()), 22);
1357 0 : EXPECT_TRUE(s.current.getEnergy() >= 0.);
1358 0 : EXPECT_TRUE(s.current.getEnergy() < Ep);
1359 :
1360 :
1361 0 : double Etot = c.current.getEnergy();
1362 0 : for (int j = 0; j < c.secondaries.size(); j++) {
1363 0 : s = *c.secondaries[j];
1364 0 : Etot += s.current.getEnergy();
1365 : }
1366 0 : EXPECT_NEAR(Ep, Etot, 1e-9);
1367 420 : }
1368 : }
1369 2 : }
1370 :
1371 2 : TEST(EMInverseComptonScattering, interactionTag) {
1372 2 : EMInverseComptonScattering m(new CMB());
1373 :
1374 : // test default interactionTag
1375 1 : EXPECT_TRUE(m.getInteractionTag() == "EMIC");
1376 :
1377 : // test secondary tag
1378 1 : m.setHavePhotons(true);
1379 2 : Candidate c(11, 1 * PeV);
1380 1 : m.performInteraction(&c);
1381 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMIC");
1382 :
1383 : // test custom tag
1384 1 : m.setInteractionTag("myTag");
1385 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
1386 1 : }
1387 :
1388 : // SynchrotronRadiation -------------------------------------------------
1389 1 : TEST(SynchrotronRadiation, interactionTag) {
1390 1 : SynchrotronRadiation s(1 * muG, true);
1391 :
1392 : // test default interactionTag
1393 1 : EXPECT_TRUE(s.getInteractionTag() == "SYN");
1394 :
1395 : // test secondary tag
1396 2 : Candidate c(11, 10 * PeV);
1397 1 : c.setCurrentStep(1 * pc);
1398 1 : s.process(&c);
1399 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "SYN");
1400 :
1401 : // test custom tag
1402 1 : s.setInteractionTag("myTag");
1403 1 : EXPECT_TRUE(s.getInteractionTag() == "myTag");
1404 1 : }
1405 :
1406 1 : TEST(SynchrotronRadiation, simpleTestRMS) {
1407 : // test initialisation with B_rms
1408 :
1409 : // check default values
1410 1 : SynchrotronRadiation sync;
1411 :
1412 1 : EXPECT_EQ(sync.getBrms(), 0);
1413 1 : EXPECT_FALSE(sync.getHavePhotons());
1414 1 : EXPECT_EQ(sync.getThinning(), 0);
1415 1 : EXPECT_EQ(sync.getLimit(), 0.1);
1416 1 : EXPECT_EQ(sync.getMaximumSamples(), 0);
1417 1 : EXPECT_EQ(sync.getSecondaryThreshold(), 1 * MeV);
1418 :
1419 : // init with custom values
1420 1 : double b = 1 * muG;
1421 1 : double thinning = 0.23;
1422 1 : int samples = 4;
1423 1 : double limit = 0.123;
1424 2 : SynchrotronRadiation sync2(b, true, thinning, samples, limit);
1425 :
1426 1 : EXPECT_EQ(sync2.getBrms(), b);
1427 1 : EXPECT_TRUE(sync2.getHavePhotons());
1428 1 : EXPECT_EQ(sync2.getThinning(), thinning);
1429 1 : EXPECT_EQ(sync2.getLimit(), limit);
1430 1 : EXPECT_EQ(sync2.getMaximumSamples(), samples);
1431 1 : EXPECT_EQ(sync2.getSecondaryThreshold(), 1 * MeV);
1432 1 : }
1433 :
1434 1 : TEST(SynchrotronRadiation, simpleTestField) {
1435 : // test initialisation with field
1436 :
1437 : // check default values
1438 : Vector3d b(0, 0, 1 * muG);
1439 1 : ref_ptr<MagneticField> field = new UniformMagneticField(b);
1440 1 : SynchrotronRadiation sync(field);
1441 :
1442 1 : EXPECT_EQ(sync.getBrms(), 0);
1443 1 : EXPECT_FALSE(sync.getHavePhotons());
1444 1 : EXPECT_EQ(sync.getThinning(), 0);
1445 1 : EXPECT_EQ(sync.getLimit(), 0.1);
1446 1 : EXPECT_EQ(sync.getMaximumSamples(), 0);
1447 1 : EXPECT_EQ(sync.getSecondaryThreshold(), 1 * MeV);
1448 1 : Vector3d fieldAtPosition = sync.getField() -> getField(Vector3d(1, 2 , 3));
1449 1 : EXPECT_EQ(fieldAtPosition.getR(), b.getR());
1450 :
1451 : // init with custom values
1452 1 : double thinning = 0.23;
1453 1 : int samples = 4;
1454 1 : double limit = 0.123;
1455 2 : SynchrotronRadiation sync2(field, true, thinning, samples, limit);
1456 :
1457 1 : EXPECT_EQ(sync2.getBrms(), 0);
1458 1 : EXPECT_TRUE(sync2.getHavePhotons());
1459 1 : EXPECT_EQ(sync2.getThinning(), thinning);
1460 1 : EXPECT_EQ(sync2.getLimit(), limit);
1461 1 : EXPECT_EQ(sync2.getMaximumSamples(), samples);
1462 1 : EXPECT_EQ(sync2.getSecondaryThreshold(), 1 * MeV);
1463 1 : fieldAtPosition = sync2.getField() -> getField(Vector3d(1, 2 , 3));
1464 1 : EXPECT_EQ(fieldAtPosition.getR(), b.getR());
1465 2 : }
1466 :
1467 1 : TEST(SynchrotronRadiation, getSetFunctions) {
1468 1 : SynchrotronRadiation sync;
1469 :
1470 : // have photons
1471 1 : sync.setHavePhotons(true);
1472 1 : EXPECT_TRUE(sync.getHavePhotons());
1473 :
1474 : // Brms
1475 1 : sync.setBrms(5 * muG);
1476 1 : EXPECT_EQ(sync.getBrms(), 5 * muG);
1477 :
1478 : // thinning
1479 1 : sync.setThinning(0.345);
1480 1 : EXPECT_EQ(sync.getThinning(), 0.345);
1481 :
1482 : // limit
1483 1 : sync.setLimit(0.234);
1484 1 : EXPECT_EQ(sync.getLimit(), 0.234);
1485 :
1486 : // max samples
1487 1 : sync.setMaximumSamples(12345);
1488 1 : EXPECT_EQ(sync.getMaximumSamples(), 12345);
1489 :
1490 : // field
1491 : Vector3d b(1,2,3);
1492 1 : ref_ptr<MagneticField> field = new UniformMagneticField(b);
1493 1 : sync.setField(field);
1494 2 : EXPECT_TRUE(field == sync.getField()); // same pointer
1495 :
1496 : // secondary threshold
1497 1 : sync.setSecondaryThreshold(1 * eV);
1498 1 : EXPECT_EQ(sync.getSecondaryThreshold(), 1 * eV);
1499 1 : }
1500 :
1501 1 : TEST(SynchrotronRadiation, energyLoss) {
1502 : double brms = 1 * muG;
1503 : double step = 1 * kpc;
1504 1 : SynchrotronRadiation sync(brms, false);
1505 :
1506 : double dE, lf, Rg, dEdx;
1507 1 : Candidate c(11);
1508 1 : c.setCurrentStep(step);
1509 1 : c.setNextStep(step);
1510 : double charge = eplus;
1511 :
1512 : // 1 GeV
1513 1 : c.current.setEnergy(1 * GeV);
1514 1 : lf = c.current.getLorentzFactor();
1515 : Rg = 1 * GeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1516 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1517 1 : dE = dEdx * step;
1518 1 : sync.process(&c);
1519 1 : EXPECT_NEAR(1 * GeV - c.current.getEnergy(), dE, 0.01 * dE);
1520 :
1521 : // 100 GeV
1522 1 : c.current.setEnergy(100 * GeV);
1523 1 : lf = c.current.getLorentzFactor();
1524 : Rg = 100 * GeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1525 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1526 1 : dE = dEdx * step;
1527 1 : sync.process(&c);
1528 1 : EXPECT_NEAR(100 * GeV - c.current.getEnergy(), dE, 0.01 * dE);
1529 :
1530 : // 10 TeV
1531 1 : c.current.setEnergy(10 * TeV);
1532 1 : lf = c.current.getLorentzFactor();
1533 : Rg = 10 * TeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1534 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1535 1 : dE = dEdx * step;
1536 1 : sync.process(&c);
1537 1 : EXPECT_NEAR(10 * TeV - c.current.getEnergy(), dE, 0.01 * dE);
1538 :
1539 : // 1 PeV
1540 1 : c.current.setEnergy(1 * PeV);
1541 1 : lf = c.current.getLorentzFactor();
1542 : Rg = 1 * PeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1543 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1544 1 : dE = dEdx * step;
1545 1 : sync.process(&c);
1546 1 : EXPECT_NEAR(1 * PeV - c.current.getEnergy(), dE, 0.01 * dE);
1547 1 : }
1548 :
1549 1 : TEST(SynchrotronRadiation, PhotonEnergy) {
1550 : double brms = 1 * muG;
1551 1 : SynchrotronRadiation sync(brms, true);
1552 1 : sync.setSecondaryThreshold(0.); // allow all secondaries for testing
1553 1 : sync.setMaximumSamples(1000); // reduce the amount of generated secondaries
1554 :
1555 : double E = 1 * TeV;
1556 1 : Candidate c(11, E);
1557 1 : c.setCurrentStep(10 * pc);
1558 1 : c.setNextStep(10 * pc);
1559 :
1560 1 : double lf = c.current.getLorentzFactor();
1561 : double Rg = E / eplus / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1562 1 : double Ecrit = 3. / 4 * h_planck / M_PI * c_light * pow(lf, 3) / Rg;
1563 :
1564 1 : sync.process(&c);
1565 1 : EXPECT_TRUE(c.secondaries.size() > 0); // must have secondaries
1566 :
1567 : // check avg energy of the secondary photons
1568 : double Esec = 0;
1569 : double weightSum = 0;
1570 1001 : for (size_t i = 0; i < c.secondaries.size(); i++) {
1571 1000 : double weight = c.secondaries[i]->getWeight();
1572 1000 : Esec += c.secondaries[i] -> current.getEnergy()*weight;
1573 1000 : weightSum += weight;
1574 : }
1575 1 : Esec /= weightSum;
1576 :
1577 1 : EXPECT_NEAR(Esec, Ecrit, Ecrit);
1578 1 : }
1579 :
1580 0 : int main(int argc, char **argv) {
1581 0 : ::testing::InitGoogleTest(&argc, argv);
1582 0 : return RUN_ALL_TESTS();
1583 : }
1584 :
1585 : } // namespace crpropa
|
