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 :
18 : #include <fstream>
19 :
20 : namespace crpropa {
21 :
22 : // ElectronPairProduction -----------------------------------------------------
23 1 : TEST(ElectronPairProduction, allBackgrounds) {
24 : // Test if interaction data files are loaded.
25 1 : ref_ptr<PhotonField> cmb = new CMB();
26 1 : ElectronPairProduction epp(cmb);
27 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
28 1 : epp.setPhotonField(irb);
29 1 : irb = new IRB_Stecker05();
30 1 : epp.setPhotonField(irb);
31 1 : irb = new IRB_Franceschini08();
32 1 : epp.setPhotonField(irb);
33 1 : irb = new IRB_Finke10();
34 1 : epp.setPhotonField(irb);
35 1 : irb = new IRB_Dominguez11();
36 1 : epp.setPhotonField(irb);
37 1 : irb = new IRB_Gilmore12();
38 1 : epp.setPhotonField(irb);
39 1 : irb = new IRB_Stecker16_upper();
40 1 : epp.setPhotonField(irb);
41 1 : irb = new IRB_Stecker16_lower();
42 1 : epp.setPhotonField(irb);
43 1 : irb = new IRB_Finke22();
44 2 : epp.setPhotonField(irb);
45 2 : }
46 :
47 1 : TEST(ElectronPairProduction, energyDecreasing) {
48 : // Test if energy loss occurs for protons with energies from 1e15 - 1e23 eV.
49 1 : Candidate c;
50 1 : c.setCurrentStep(2 * Mpc);
51 1 : c.current.setId(nucleusId(1, 1)); // proton
52 :
53 1 : ref_ptr<PhotonField> cmb = new CMB();
54 1 : ElectronPairProduction epp1(cmb);
55 81 : for (int i = 0; i < 80; i++) {
56 80 : double E = pow(10, 15 + i * 0.1) * eV;
57 80 : c.current.setEnergy(E);
58 80 : epp1.process(&c);
59 80 : EXPECT_LE(c.current.getEnergy(), E);
60 : }
61 :
62 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
63 2 : ElectronPairProduction epp2(irb);
64 81 : for (int i = 0; i < 80; i++) {
65 80 : double E = pow(10, 15 + i * 0.1) * eV;
66 80 : c.current.setEnergy(E);
67 80 : epp2.process(&c);
68 80 : EXPECT_LE(c.current.getEnergy(), E);
69 : }
70 2 : }
71 :
72 1 : TEST(ElectronPairProduction, belowEnergyTreshold) {
73 : // Test if nothing happens below 1e15 eV.
74 1 : ref_ptr<PhotonField> cmb = new CMB();
75 1 : ElectronPairProduction epp(cmb);
76 1 : Candidate c(nucleusId(1, 1), 1E14 * eV);
77 1 : epp.process(&c);
78 1 : EXPECT_DOUBLE_EQ(1E14 * eV, c.current.getEnergy());
79 2 : }
80 :
81 1 : TEST(ElectronPairProduction, thisIsNotNucleonic) {
82 : // Test if non-nuclei are skipped.
83 1 : ref_ptr<PhotonField> cmb = new CMB();
84 1 : ElectronPairProduction epp(cmb);
85 1 : Candidate c(11, 1E20 * eV); // electron
86 1 : epp.process(&c);
87 1 : EXPECT_DOUBLE_EQ(1E20 * eV, c.current.getEnergy());
88 2 : }
89 :
90 2 : TEST(ElectronPairProduction, valuesCMB) {
91 : // Test if energy loss corresponds to the data table.
92 : std::vector<double> x;
93 : std::vector<double> y;
94 1 : std::ifstream infile(getDataPath("pair_CMB.txt").c_str());
95 1 : while (infile.good()) {
96 0 : if (infile.peek() != '#') {
97 : double a, b;
98 : infile >> a >> b;
99 0 : if (infile) {
100 0 : x.push_back(a * eV);
101 0 : y.push_back(b * eV / Mpc);
102 : }
103 : }
104 0 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
105 : }
106 1 : infile.close();
107 :
108 1 : Candidate c;
109 1 : c.setCurrentStep(1 * Mpc);
110 1 : c.current.setId(nucleusId(1, 1)); // proton
111 1 : ref_ptr<PhotonField> cmb = new CMB();
112 :
113 1 : ElectronPairProduction epp(cmb);
114 1 : for (int i = 0; i < x.size(); i++) {
115 0 : c.current.setEnergy(x[i]);
116 0 : epp.process(&c);
117 0 : double dE = x[i] - c.current.getEnergy();
118 0 : double dE_table = y[i] * 1 * Mpc;
119 0 : EXPECT_NEAR(dE_table, dE, 1e-12);
120 : }
121 2 : }
122 :
123 1 : TEST(ElectronPairProduction, interactionTag) {
124 :
125 1 : ref_ptr<PhotonField> cmb = new CMB();
126 0 : ElectronPairProduction epp(cmb);
127 :
128 : // test the default interaction tag
129 1 : EXPECT_TRUE(epp.getInteractionTag() == "EPP");
130 :
131 : // test changing the interaction tag
132 1 : epp.setInteractionTag("myTag");
133 1 : EXPECT_TRUE(epp.getInteractionTag() == "myTag");
134 :
135 : // test the tag of produced secondaries
136 2 : Candidate c;
137 1 : c.setCurrentStep(1 * Gpc);
138 1 : c.current.setId(nucleusId(1,1));
139 1 : c.current.setEnergy(100 * EeV);
140 1 : epp.setHaveElectrons(true);
141 1 : epp.process(&c);
142 :
143 1 : std::string secondaryTag = c.secondaries[0] -> getTagOrigin();
144 1 : EXPECT_TRUE(secondaryTag == "myTag");
145 2 : }
146 :
147 2 : TEST(ElectronPairProduction, valuesIRB) {
148 : // Test if energy loss corresponds to the data table.
149 : std::vector<double> x;
150 : std::vector<double> y;
151 1 : std::ifstream infile(getDataPath("pairIRB.txt").c_str());
152 1 : while (infile.good()) {
153 0 : if (infile.peek() != '#') {
154 : double a, b;
155 : infile >> a >> b;
156 0 : if (infile) {
157 0 : x.push_back(a * eV);
158 0 : y.push_back(b * eV / Mpc);
159 : }
160 : }
161 0 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
162 : }
163 1 : infile.close();
164 :
165 1 : Candidate c;
166 1 : c.setCurrentStep(1 * Mpc);
167 1 : c.current.setId(nucleusId(1, 1)); // proton
168 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
169 :
170 1 : ElectronPairProduction epp(irb);
171 1 : for (int i = 0; i < x.size(); i++) {
172 0 : c.current.setEnergy(x[i]);
173 0 : epp.process(&c);
174 0 : double dE = x[i] - c.current.getEnergy();
175 0 : double dE_table = y[i] * 1 * Mpc;
176 0 : EXPECT_NEAR(dE, dE_table, 1e-12);
177 : }
178 2 : }
179 :
180 : // NuclearDecay ---------------------------------------------------------------
181 1 : TEST(NuclearDecay, scandium44) {
182 : // Test beta+ decay of 44Sc to 44Ca.
183 : // This test can stochastically fail.
184 1 : NuclearDecay d(true, true);
185 1 : Candidate c(nucleusId(44, 21), 1E18 * eV);
186 1 : c.setCurrentStep(100 * Mpc);
187 1 : double gamma = c.current.getLorentzFactor();
188 1 : d.process(&c);
189 :
190 : // expected decay product: 44Ca
191 1 : EXPECT_EQ(nucleusId(44, 20), c.current.getId());
192 :
193 : // expect Lorentz factor to be conserved
194 1 : EXPECT_DOUBLE_EQ(gamma, c.current.getLorentzFactor());
195 :
196 : // expect at least two secondaries: positron + electron neutrino
197 1 : EXPECT_GE(c.secondaries.size(), 2);
198 1 : }
199 :
200 1 : TEST(NuclearDecay, lithium4) {
201 : // Test proton dripping of Li-4 to He-3
202 : // This test can stochastically fail
203 1 : NuclearDecay d;
204 1 : Candidate c(nucleusId(4, 3), 4 * EeV);
205 1 : c.setCurrentStep(100 * Mpc);
206 1 : d.process(&c);
207 :
208 : // expected decay product: He-3
209 1 : EXPECT_EQ(nucleusId(3, 2), c.current.getId());
210 :
211 : // expected secondary: proton
212 1 : EXPECT_EQ(1, c.secondaries.size());
213 2 : Candidate c1 = *c.secondaries[0];
214 1 : EXPECT_EQ(nucleusId(1, 1), c1.current.getId());
215 1 : EXPECT_EQ(1 * EeV, c1.current.getEnergy());
216 1 : }
217 :
218 1 : TEST(NuclearDecay, helium5) {
219 : // Test neutron dripping of He-5 to He-4.
220 : // This test can stochastically fail.
221 1 : NuclearDecay d;
222 1 : Candidate c(nucleusId(5, 2), 5 * EeV);
223 1 : c.setCurrentStep(100 * Mpc);
224 1 : d.process(&c);
225 :
226 : // expected primary: He-4
227 1 : EXPECT_EQ(nucleusId(4, 2), c.current.getId());
228 1 : EXPECT_EQ(4, c.current.getEnergy() / EeV);
229 :
230 : // expected secondary: neutron
231 2 : Candidate c2 = *c.secondaries[0];
232 1 : EXPECT_EQ(nucleusId(1, 0), c2.current.getId());
233 1 : EXPECT_EQ(1, c2.current.getEnergy() / EeV);
234 1 : }
235 :
236 1 : TEST(NuclearDecay, limitNextStep) {
237 : // Test if next step is limited in case of a neutron.
238 1 : NuclearDecay decay;
239 1 : Candidate c(nucleusId(1, 0), 10 * EeV);
240 1 : c.setNextStep(std::numeric_limits<double>::max());
241 1 : decay.process(&c);
242 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
243 1 : }
244 :
245 1 : TEST(NuclearDecay, allChannelsWorking) {
246 : // Test if all nuclear decays are working.
247 1 : NuclearDecay d;
248 1 : Candidate c;
249 :
250 1 : std::ifstream infile(getDataPath("nuclear_decay.txt").c_str());
251 951 : while (infile.good()) {
252 950 : if (infile.peek() != '#') {
253 : int Z, N, channel, foo;
254 948 : infile >> Z >> N >> channel >> foo;
255 948 : c.current.setId(nucleusId(Z + N, Z));
256 948 : c.current.setEnergy(80 * EeV);
257 948 : d.performInteraction(&c, channel);
258 : }
259 950 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
260 : }
261 1 : infile.close();
262 1 : }
263 :
264 1 : TEST(NuclearDecay, secondaries) {
265 : // Test if all types of secondaries are produced.
266 1 : NuclearDecay d;
267 1 : d.setHaveElectrons(true);
268 1 : d.setHaveNeutrinos(true);
269 1 : d.setHavePhotons(true);
270 1 : Candidate c;
271 :
272 : // He-8 --> Li-8 + e- + neutrino
273 : // additional photon emitted with 84% probability
274 : // --> expect at least 1 photon out of 10 decays
275 11 : for (int i = 0; i < 10; ++i) {
276 10 : c.current.setId(nucleusId(8, 2));
277 10 : c.current.setEnergy(5 * EeV);
278 10 : d.performInteraction(&c, 10000);
279 : }
280 :
281 : // count number of secondaries
282 1 : size_t nElectrons = 0;
283 1 : size_t nNeutrinos = 0;
284 1 : size_t nPhotons = 0;
285 :
286 28 : for(size_t i = 0; i < c.secondaries.size(); ++i) {
287 27 : int id = (*c.secondaries[i]).current.getId();
288 27 : if (id == 22) nPhotons++;
289 27 : if (id == 11) nElectrons++;
290 27 : if (id == -12) nNeutrinos++;
291 : }
292 :
293 1 : EXPECT_EQ(nElectrons, 10);
294 1 : EXPECT_EQ(nNeutrinos, 10);
295 1 : EXPECT_GE(nPhotons, 1);
296 1 : }
297 :
298 1 : TEST(NuclearDecay, thisIsNotNucleonic) {
299 : // Test if nothing happens to an electron
300 1 : NuclearDecay decay;
301 1 : Candidate c(11, 10 * EeV);
302 1 : c.setNextStep(std::numeric_limits<double>::max());
303 1 : decay.process(&c);
304 1 : EXPECT_EQ(11, c.current.getId());
305 1 : EXPECT_EQ(10 * EeV, c.current.getEnergy());
306 1 : }
307 :
308 1 : TEST(NuclearDecay, interactionTag) {
309 : // test default interaction tag
310 1 : NuclearDecay decay;
311 1 : EXPECT_TRUE(decay.getInteractionTag() == "ND");
312 :
313 : // test secondary tag
314 1 : decay.setHaveElectrons(true);
315 2 : Candidate c(nucleusId(8,2), 5 * EeV);
316 1 : decay.performInteraction(&c, 10000);
317 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "ND");
318 :
319 : // test custom tags
320 1 : decay.setInteractionTag("myTag");
321 1 : EXPECT_TRUE(decay.getInteractionTag() == "myTag");
322 1 : }
323 :
324 : // PhotoDisintegration --------------------------------------------------------
325 1 : TEST(PhotoDisintegration, allBackgrounds) {
326 : // Test if interaction data files are loaded.
327 1 : ref_ptr<PhotonField> cmb = new CMB();
328 1 : PhotoDisintegration pd(cmb);
329 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
330 1 : pd.setPhotonField(irb);
331 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
332 1 : pd.setPhotonField(urb);
333 1 : irb = new IRB_Stecker05();
334 1 : pd.setPhotonField(irb);
335 1 : irb = new IRB_Franceschini08();
336 1 : pd.setPhotonField(irb);
337 1 : irb = new IRB_Finke10();
338 1 : pd.setPhotonField(irb);
339 1 : irb = new IRB_Dominguez11();
340 1 : pd.setPhotonField(irb);
341 1 : irb = new IRB_Gilmore12();
342 1 : pd.setPhotonField(irb);
343 1 : irb = new IRB_Stecker16_upper();
344 1 : pd.setPhotonField(irb);
345 1 : irb = new IRB_Stecker16_lower();
346 1 : pd.setPhotonField(irb);
347 1 : irb = new IRB_Finke22();
348 1 : pd.setPhotonField(irb);
349 1 : urb = new URB_Nitu21();
350 2 : pd.setPhotonField(urb);
351 2 : }
352 :
353 1 : TEST(PhotoDisintegration, carbon) {
354 : // Test if a 100 EeV C-12 nucleus photo-disintegrates (at least once) over a distance of 1 Gpc.
355 : // This test can stochastically fail.
356 1 : ref_ptr<PhotonField> cmb = new CMB();
357 1 : PhotoDisintegration pd(cmb);
358 1 : Candidate c;
359 1 : int id = nucleusId(12, 6);
360 1 : c.current.setId(id);
361 1 : c.current.setEnergy(100 * EeV);
362 1 : c.setCurrentStep(1000 * Mpc);
363 1 : pd.process(&c);
364 :
365 1 : EXPECT_TRUE(c.current.getEnergy() < 100 * EeV);
366 : // energy loss
367 1 : EXPECT_TRUE(c.secondaries.size() > 0);
368 : // secondaries produced
369 :
370 1 : double E = c.current.getEnergy();
371 1 : id = c.current.getId();
372 1 : int A = massNumber(id);
373 1 : int Z = chargeNumber(id);
374 :
375 4 : for (int i = 0; i < c.secondaries.size(); i++) {
376 3 : E += (*c.secondaries[i]).current.getEnergy();
377 3 : id = (*c.secondaries[i]).current.getId();
378 3 : A += massNumber(id);
379 3 : Z += chargeNumber(id);
380 : }
381 1 : EXPECT_EQ(12, A);
382 : // nucleon number conserved
383 1 : EXPECT_EQ(6, Z);
384 : // proton number conserved
385 1 : EXPECT_DOUBLE_EQ(100 * EeV, E);
386 : // energy conserved
387 2 : }
388 :
389 1 : TEST(PhotoDisintegration, iron) {
390 : // Test if a 200 EeV Fe-56 nucleus photo-disintegrates (at least once) over a distance of 1 Gpc.
391 : // This test can stochastically fail.
392 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
393 1 : PhotoDisintegration pd(irb);
394 1 : Candidate c;
395 1 : int id = nucleusId(56, 26);
396 1 : c.current.setId(id);
397 1 : c.current.setEnergy(200 * EeV);
398 1 : c.setCurrentStep(1000 * Mpc);
399 1 : pd.process(&c);
400 :
401 : // expect energy loss
402 1 : EXPECT_TRUE(c.current.getEnergy() < 200 * EeV);
403 :
404 : // expect secondaries produced
405 1 : EXPECT_TRUE(c.secondaries.size() > 0);
406 :
407 1 : double E = c.current.getEnergy();
408 1 : id = c.current.getId();
409 1 : int A = massNumber(id);
410 1 : int Z = chargeNumber(id);
411 :
412 43 : for (int i = 0; i < c.secondaries.size(); i++) {
413 42 : E += (*c.secondaries[i]).current.getEnergy();
414 42 : id = (*c.secondaries[i]).current.getId();
415 42 : A += massNumber(id);
416 42 : Z += chargeNumber(id);
417 : }
418 :
419 : // nucleon number conserved
420 1 : EXPECT_EQ(56, A);
421 :
422 : // proton number conserved (no decay active)
423 1 : EXPECT_EQ(26, Z);
424 :
425 : // energy conserved
426 1 : EXPECT_DOUBLE_EQ(200 * EeV, E);
427 2 : }
428 :
429 1 : TEST(PhotoDisintegration, thisIsNotNucleonic) {
430 : // Test that nothing happens to an electron.
431 1 : ref_ptr<PhotonField> cmb = new CMB();
432 1 : PhotoDisintegration pd(cmb);
433 1 : Candidate c;
434 1 : c.setCurrentStep(1 * Mpc);
435 1 : c.current.setId(11); // electron
436 1 : c.current.setEnergy(10 * EeV);
437 1 : pd.process(&c);
438 1 : EXPECT_EQ(11, c.current.getId());
439 1 : EXPECT_EQ(10 * EeV, c.current.getEnergy());
440 2 : }
441 :
442 1 : TEST(PhotoDisintegration, limitNextStep) {
443 : // Test if the interaction limits the next propagation step.
444 1 : ref_ptr<PhotonField> cmb = new CMB();
445 1 : PhotoDisintegration pd(cmb);
446 1 : Candidate c;
447 1 : c.setNextStep(std::numeric_limits<double>::max());
448 1 : c.current.setId(nucleusId(4, 2));
449 1 : c.current.setEnergy(200 * EeV);
450 1 : pd.process(&c);
451 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
452 2 : }
453 :
454 1 : TEST(PhotoDisintegration, allIsotopes) {
455 : // Test if all isotopes are handled.
456 1 : ref_ptr<PhotonField> cmb = new CMB();
457 1 : PhotoDisintegration pd1(cmb);
458 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
459 1 : PhotoDisintegration pd2(irb);
460 1 : Candidate c;
461 1 : c.setCurrentStep(10 * Mpc);
462 :
463 27 : for (int Z = 1; Z <= 26; Z++) {
464 806 : for (int N = 1; N <= 30; N++) {
465 :
466 780 : c.current.setId(nucleusId(Z + N, Z));
467 780 : c.current.setEnergy(80 * EeV);
468 780 : pd1.process(&c);
469 :
470 780 : c.current.setId(nucleusId(Z + N, Z));
471 780 : c.current.setEnergy(80 * EeV);
472 780 : pd2.process(&c);
473 : }
474 : }
475 3 : }
476 :
477 1 : TEST(Photodisintegration, updateParticleParentProperties) { // Issue: #204
478 1 : ref_ptr<PhotonField> cmb = new CMB();
479 1 : PhotoDisintegration pd(cmb);
480 :
481 1 : Candidate c(nucleusId(56,26), 500 * EeV, Vector3d(1 * Mpc, 0, 0));
482 :
483 1 : pd.performInteraction(&c, 1);
484 : // the candidates parent is the original particle
485 1 : EXPECT_EQ(c.created.getId(), nucleusId(56,26));
486 :
487 1 : pd.performInteraction(&c, 1);
488 : // now it has to be changed
489 1 : EXPECT_NE(c.created.getId(), nucleusId(56,26));
490 2 : }
491 :
492 1 : TEST(PhotoDisintegration, interactionTag) {
493 1 : PhotoDisintegration pd(new CMB());
494 :
495 : // test default interactionTag
496 1 : EXPECT_TRUE(pd.getInteractionTag() == "PD");
497 :
498 : // test secondary tag
499 1 : pd.setHavePhotons(true);
500 2 : Candidate c(nucleusId(56,26), 500 * EeV);
501 1 : c.setCurrentStep(1 * Gpc);
502 1 : pd.process(&c);
503 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "PD");
504 :
505 : // test custom tag
506 1 : pd.setInteractionTag("myTag");
507 1 : EXPECT_TRUE(pd.getInteractionTag() == "myTag");
508 1 : }
509 :
510 : // ElasticScattering ----------------------------------------------------------
511 1 : TEST(ElasticScattering, allBackgrounds) {
512 : // Test if interaction data files are loaded.
513 1 : ref_ptr<PhotonField> cmb = new CMB();
514 1 : ElasticScattering scattering(cmb);
515 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
516 1 : scattering.setPhotonField(irb);
517 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
518 2 : scattering.setPhotonField(urb);
519 2 : }
520 :
521 1 : TEST(ElasticScattering, secondaries) {
522 : // Test the creation of cosmic ray photons.
523 : // This test can stochastically fail.
524 1 : ref_ptr<PhotonField> cmb = new CMB();
525 1 : ElasticScattering scattering(cmb);
526 1 : Candidate c;
527 1 : int id = nucleusId(12, 6);
528 1 : c.current.setId(id);
529 1 : c.current.setEnergy(200 * EeV);
530 1 : c.setCurrentStep(400 * Mpc);
531 1 : scattering.process(&c);
532 :
533 1 : EXPECT_GT(c.secondaries.size(), 0);
534 :
535 7 : for (int i = 0; i < c.secondaries.size(); i++) {
536 6 : int id = (*c.secondaries[i]).current.getId();
537 6 : EXPECT_EQ(id, 22);
538 6 : double energy = (*c.secondaries[i]).current.getEnergy();
539 6 : EXPECT_GT(energy, 0);
540 6 : EXPECT_LT(energy, 200 * EeV);
541 : }
542 2 : }
543 :
544 : // PhotoPionProduction --------------------------------------------------------
545 1 : TEST(PhotoPionProduction, allBackgrounds) {
546 : // Test if all interaction data files can be loaded.
547 1 : ref_ptr<PhotonField> cmb = new CMB();
548 1 : PhotoPionProduction ppp(cmb);
549 1 : ref_ptr<PhotonField> irb = new IRB_Kneiske04();
550 1 : ppp.setPhotonField(irb);
551 1 : irb = new IRB_Stecker05();
552 1 : ppp.setPhotonField(irb);
553 1 : irb = new IRB_Franceschini08();
554 1 : ppp.setPhotonField(irb);
555 1 : irb = new IRB_Finke10();
556 1 : ppp.setPhotonField(irb);
557 1 : irb = new IRB_Dominguez11();
558 1 : ppp.setPhotonField(irb);
559 1 : irb = new IRB_Gilmore12();
560 1 : ppp.setPhotonField(irb);
561 1 : irb = new IRB_Stecker16_upper();
562 1 : ppp.setPhotonField(irb);
563 1 : irb = new IRB_Stecker16_lower();
564 1 : ppp.setPhotonField(irb);
565 1 : irb = new IRB_Finke22();
566 1 : ppp.setPhotonField(irb);
567 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
568 1 : ppp.setPhotonField(urb);
569 1 : urb = new URB_Nitu21();
570 2 : ppp.setPhotonField(urb);
571 2 : }
572 :
573 1 : TEST(PhotoPionProduction, proton) {
574 : // Test photopion interaction for 100 EeV proton.
575 : // This test can stochastically fail.
576 1 : ref_ptr<PhotonField> cmb = new CMB();
577 1 : PhotoPionProduction ppp(cmb);
578 1 : Candidate c(nucleusId(1, 1), 100 * EeV);
579 1 : c.setCurrentStep(1000 * Mpc);
580 1 : ppp.process(&c);
581 :
582 : // expect energy loss
583 1 : EXPECT_LT(c.current.getEnergy(), 100. * EeV);
584 :
585 : // expect nucleon number conservation
586 1 : EXPECT_EQ(1, massNumber(c.current.getId()));
587 :
588 : // expect no (nucleonic) secondaries
589 1 : EXPECT_EQ(0, c.secondaries.size());
590 2 : }
591 :
592 1 : TEST(PhotoPionProduction, helium) {
593 : // Test photo-pion interaction for 400 EeV He nucleus.
594 : // This test can stochastically fail.
595 1 : ref_ptr<PhotonField> cmb = new CMB();
596 1 : PhotoPionProduction ppp(cmb);
597 1 : Candidate c;
598 1 : c.current.setId(nucleusId(4, 2));
599 1 : c.current.setEnergy(400. * EeV);
600 1 : c.setCurrentStep(1000 * Mpc);
601 1 : ppp.process(&c);
602 1 : EXPECT_LT(c.current.getEnergy(), 400. * EeV);
603 1 : int id = c.current.getId();
604 1 : EXPECT_TRUE(massNumber(id) < 4);
605 1 : EXPECT_TRUE(c.secondaries.size() > 0);
606 2 : }
607 :
608 1 : TEST(PhotoPionProduction, thisIsNotNucleonic) {
609 : // Test if nothing happens to an electron.
610 1 : ref_ptr<PhotonField> cmb = new CMB();
611 1 : PhotoPionProduction ppp(cmb);
612 1 : Candidate c;
613 1 : c.current.setId(11); // electron
614 1 : c.current.setEnergy(10 * EeV);
615 1 : c.setCurrentStep(100 * Mpc);
616 1 : ppp.process(&c);
617 1 : EXPECT_EQ(11, c.current.getId());
618 1 : EXPECT_EQ(10 * EeV, c.current.getEnergy());
619 2 : }
620 :
621 1 : TEST(PhotoPionProduction, limitNextStep) {
622 : // Test if the interaction limits the next propagation step.
623 1 : ref_ptr<PhotonField> cmb = new CMB();
624 1 : PhotoPionProduction ppp(cmb);
625 1 : Candidate c(nucleusId(1, 1), 200 * EeV);
626 1 : c.setNextStep(std::numeric_limits<double>::max());
627 1 : ppp.process(&c);
628 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
629 2 : }
630 :
631 1 : TEST(PhotoPionProduction, secondaries) {
632 : // Test photo-pion interaction for 100 EeV proton.
633 : // This test can stochastically fail.
634 1 : ref_ptr<PhotonField> cmb = new CMB();
635 1 : PhotoPionProduction ppp(cmb, true, true, true);
636 1 : Candidate c(nucleusId(1, 1), 100 * EeV);
637 1 : c.setCurrentStep(1000 * Mpc);
638 1 : ppp.process(&c);
639 : // there should be secondaries
640 1 : EXPECT_GT(c.secondaries.size(), 1);
641 2 : }
642 :
643 1 : TEST(PhotoPionProduction, sampling) {
644 : // Specific test of photon sampling of photo-pion production
645 : // by testing the calculated pEpsMax for CMB(), also indirectly
646 : // testing epsMinInteraction and logSampling (default).
647 1 : ref_ptr<PhotonField> cmb = new CMB(); //create CMB instance
648 : double energy = 1.e10; //1e10 GeV
649 : bool onProton = true; //proton
650 : double z = 0; //no redshift
651 1 : PhotoPionProduction ppp(cmb, true, true, true);
652 1 : double correctionFactor = ppp.getCorrectionFactor(); //get current correctionFactor
653 1 : double epsMin = std::max(cmb -> getMinimumPhotonEnergy(z) / eV, 0.00710614); // 0.00710614 = epsMinInteraction(onProton,energy)
654 1 : double epsMax = cmb -> getMaximumPhotonEnergy(z) / eV;
655 1 : double pEpsMax = ppp.probEpsMax(onProton, energy, z, epsMin, epsMax) / correctionFactor;
656 1 : EXPECT_DOUBLE_EQ(pEpsMax,132673934934.922);
657 2 : }
658 :
659 1 : TEST(PhotoPionProduction, interactionTag) {
660 1 : PhotoPionProduction ppp(new CMB());
661 :
662 : // test default interactionTag
663 1 : EXPECT_TRUE(ppp.getInteractionTag() == "PPP");
664 :
665 : // test secondary tag
666 1 : ppp.setHavePhotons(true);
667 2 : Candidate c(nucleusId(1,1), 100 * EeV);
668 11 : for(int i = 0; i <10; i++)
669 10 : ppp.performInteraction(&c, true);
670 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "PPP");
671 :
672 : // test custom interactionTag
673 1 : ppp.setInteractionTag("myTag");
674 1 : EXPECT_TRUE(ppp.getInteractionTag() == "myTag");
675 1 : }
676 :
677 : // Redshift -------------------------------------------------------------------
678 2 : TEST(Redshift, simpleTest) {
679 : // Test if redshift is decreased and adiabatic energy loss is applied.
680 : Redshift redshift;
681 :
682 1 : Candidate c;
683 1 : c.setRedshift(0.024);
684 1 : c.current.setEnergy(100 * EeV);
685 1 : c.setCurrentStep(1 * Mpc);
686 :
687 1 : redshift.process(&c);
688 1 : EXPECT_GT(0.024, c.getRedshift()); // expect redshift decrease
689 1 : EXPECT_GT(100, c.current.getEnergy() / EeV); // expect energy loss
690 2 : }
691 :
692 2 : TEST(Redshift, limitRedshiftDecrease) {
693 : // Test if the redshift decrease is limited to z_updated >= 0.
694 : Redshift redshift;
695 :
696 1 : Candidate c;
697 1 : c.setRedshift(0.024); // roughly corresponds to 100 Mpc
698 1 : c.setCurrentStep(150 * Mpc);
699 :
700 1 : redshift.process(&c);
701 1 : EXPECT_DOUBLE_EQ(0, c.getRedshift());
702 2 : }
703 :
704 : // EMPairProduction -----------------------------------------------------------
705 1 : TEST(EMPairProduction, allBackgrounds) {
706 : // Test if interaction data files are loaded.
707 1 : ref_ptr<PhotonField> cmb = new CMB();
708 1 : EMPairProduction em(cmb);
709 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
710 1 : em.setPhotonField(ebl);
711 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
712 1 : em.setPhotonField(urb);
713 1 : ebl = new IRB_Stecker05();
714 1 : em.setPhotonField(ebl);
715 1 : ebl = new IRB_Franceschini08();
716 1 : em.setPhotonField(ebl);
717 1 : ebl = new IRB_Finke10();
718 1 : em.setPhotonField(ebl);
719 1 : ebl = new IRB_Dominguez11();
720 1 : em.setPhotonField(ebl);
721 1 : ebl = new IRB_Gilmore12();
722 1 : em.setPhotonField(ebl);
723 1 : ebl = new IRB_Stecker16_upper();
724 1 : em.setPhotonField(ebl);
725 1 : ebl = new IRB_Stecker16_lower();
726 1 : em.setPhotonField(ebl);
727 1 : ebl = new IRB_Finke22();
728 1 : em.setPhotonField(ebl);
729 1 : urb = new URB_Fixsen11();
730 1 : em.setPhotonField(urb);
731 1 : urb = new URB_Nitu21();
732 2 : em.setPhotonField(urb);
733 2 : }
734 :
735 1 : TEST(EMPairProduction, limitNextStep) {
736 : // Test if the interaction limits the next propagation step.
737 1 : ref_ptr<PhotonField> cmb = new CMB();
738 1 : EMPairProduction m(cmb);
739 1 : Candidate c(22, 1E17 * eV);
740 1 : c.setNextStep(std::numeric_limits<double>::max());
741 1 : m.process(&c);
742 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
743 2 : }
744 :
745 1 : TEST(EMPairProduction, secondaries) {
746 : // Test if secondaries are correctly produced.
747 1 : ref_ptr<PhotonField> cmb = new CMB();
748 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
749 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
750 1 : EMPairProduction m(cmb);
751 1 : m.setHaveElectrons(true);
752 1 : m.setThinning(0.);
753 :
754 : std::vector<ref_ptr<PhotonField>> fields;
755 1 : fields.push_back(cmb);
756 1 : fields.push_back(irb);
757 1 : fields.push_back(urb);
758 :
759 : // loop over photon backgrounds
760 4 : for (int f = 0; f < fields.size(); f++) {
761 3 : m.setPhotonField(fields[f]);
762 423 : for (int i = 0; i < 140; i++) { // loop over energies Ep = (1e10 - 1e23) eV
763 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
764 420 : Candidate c(22, Ep);
765 420 : c.setCurrentStep(1e10 * Mpc);
766 :
767 420 : m.process(&c);
768 :
769 : // pass if no interaction has ocurred (no tabulated rates)
770 420 : if (c.isActive())
771 : continue;
772 :
773 : // expect 2 secondaries
774 311 : EXPECT_EQ(c.secondaries.size(), 2);
775 :
776 : // expect electron / positron with energies 0 < E < Ephoton
777 : double Etot = 0;
778 933 : for (int j = 0; j < c.secondaries.size(); j++) {
779 622 : Candidate s = *c.secondaries[j];
780 622 : EXPECT_EQ(abs(s.current.getId()), 11);
781 622 : EXPECT_GT(s.current.getEnergy(), 0);
782 622 : EXPECT_LT(s.current.getEnergy(), Ep);
783 622 : Etot += s.current.getEnergy();
784 622 : }
785 :
786 : // test energy conservation
787 311 : EXPECT_DOUBLE_EQ(Ep, Etot);
788 420 : }
789 : }
790 2 : }
791 :
792 1 : TEST(EMPairProduction, interactionTag) {
793 1 : EMPairProduction m(new CMB());
794 :
795 : // test default interactionTag
796 1 : EXPECT_TRUE(m.getInteractionTag() == "EMPP");
797 :
798 : // test secondary tag
799 1 : m.setHaveElectrons(true);
800 2 : Candidate c(22, 1 * EeV);
801 1 : m.performInteraction(&c);
802 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMPP");
803 :
804 : // test custom tag
805 1 : m.setInteractionTag("myTag");
806 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
807 1 : }
808 :
809 : // EMDoublePairProduction -----------------------------------------------------
810 1 : TEST(EMDoublePairProduction, allBackgrounds) {
811 : // Test if interaction data files are loaded.
812 1 : ref_ptr<PhotonField> cmb = new CMB();
813 1 : EMDoublePairProduction em(cmb);
814 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
815 1 : em.setPhotonField(ebl);
816 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
817 1 : em.setPhotonField(urb);
818 1 : ebl = new IRB_Stecker05();
819 1 : em.setPhotonField(ebl);
820 1 : ebl = new IRB_Franceschini08();
821 1 : em.setPhotonField(ebl);
822 1 : ebl = new IRB_Finke10();
823 1 : em.setPhotonField(ebl);
824 1 : ebl = new IRB_Dominguez11();
825 1 : em.setPhotonField(ebl);
826 1 : ebl = new IRB_Gilmore12();
827 1 : em.setPhotonField(ebl);
828 1 : ebl = new IRB_Stecker16_upper();
829 1 : em.setPhotonField(ebl);
830 1 : ebl = new IRB_Stecker16_lower();
831 1 : em.setPhotonField(ebl);
832 1 : ebl = new IRB_Finke22();
833 1 : em.setPhotonField(ebl);
834 1 : urb = new URB_Fixsen11();
835 1 : em.setPhotonField(urb);
836 1 : urb = new URB_Nitu21();
837 2 : em.setPhotonField(urb);
838 2 : }
839 :
840 1 : TEST(EMDoublePairProduction, limitNextStep) {
841 : // Test if the interaction limits the next propagation step.
842 1 : ref_ptr<PhotonField> cmb = new CMB();
843 1 : EMDoublePairProduction m(cmb);
844 1 : Candidate c(22, 1E17 * eV);
845 1 : c.setNextStep(std::numeric_limits<double>::max());
846 1 : m.process(&c);
847 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
848 2 : }
849 :
850 1 : TEST(EMDoublePairProduction, secondaries) {
851 : // Test if secondaries are correctly produced.
852 1 : ref_ptr<PhotonField> cmb = new CMB();
853 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
854 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
855 1 : EMPairProduction m(cmb);
856 1 : m.setHaveElectrons(true);
857 1 : m.setThinning(0.);
858 :
859 : std::vector<ref_ptr<PhotonField>> fields;
860 1 : fields.push_back(cmb);
861 1 : fields.push_back(irb);
862 1 : fields.push_back(urb);
863 :
864 : // loop over photon backgrounds
865 4 : for (int f = 0; f < fields.size(); f++) {
866 3 : m.setPhotonField(fields[f]);
867 :
868 : // loop over energies Ep = (1e9 - 1e23) eV
869 423 : for (int i = 0; i < 140; i++) {
870 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
871 420 : Candidate c(22, Ep);
872 420 : c.setCurrentStep(1e4 * Mpc); // use lower value so that the test can run faster
873 420 : m.process(&c);
874 :
875 : // pass if no interaction has occured (no tabulated rates)
876 420 : if (c.isActive())
877 : continue;
878 :
879 : // expect 2 secondaries (only one pair is considered)
880 249 : EXPECT_EQ(c.secondaries.size(), 2);
881 :
882 : // expect electron / positron with energies 0 < E < Ephoton
883 : double Etot = 0;
884 747 : for (int j = 0; j < c.secondaries.size(); j++) {
885 498 : Candidate s = *c.secondaries[j];
886 498 : EXPECT_EQ(abs(s.current.getId()), 11);
887 498 : EXPECT_GT(s.current.getEnergy(), 0);
888 498 : EXPECT_LT(s.current.getEnergy(), Ep);
889 498 : Etot += s.current.getEnergy();
890 498 : }
891 :
892 : // test energy conservation
893 249 : EXPECT_NEAR(Ep, Etot, 1E-9);
894 420 : }
895 : }
896 2 : }
897 :
898 1 : TEST(EMDoublePairProduction, interactionTag) {
899 1 : EMDoublePairProduction m(new CMB());
900 :
901 : // test default interactionTag
902 1 : EXPECT_TRUE(m.getInteractionTag() == "EMDP");
903 :
904 : // test secondary tag
905 1 : m.setHaveElectrons(true);
906 2 : Candidate c(22, 1 * EeV);
907 1 : m.performInteraction(&c);
908 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMDP");
909 :
910 : // test custom tag
911 1 : m.setInteractionTag("myTag");
912 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
913 1 : }
914 :
915 : // EMTripletPairProduction ----------------------------------------------------
916 1 : TEST(EMTripletPairProduction, allBackgrounds) {
917 : // Test if interaction data files are loaded.
918 1 : ref_ptr<PhotonField> cmb = new CMB();
919 1 : EMTripletPairProduction em(cmb);
920 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
921 1 : em.setPhotonField(ebl);
922 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
923 1 : em.setPhotonField(urb);
924 1 : ebl = new IRB_Stecker05();
925 1 : em.setPhotonField(ebl);
926 1 : ebl = new IRB_Franceschini08();
927 1 : em.setPhotonField(ebl);
928 1 : ebl = new IRB_Finke10();
929 1 : em.setPhotonField(ebl);
930 1 : ebl = new IRB_Dominguez11();
931 1 : em.setPhotonField(ebl);
932 1 : ebl = new IRB_Gilmore12();
933 1 : em.setPhotonField(ebl);
934 1 : ebl = new IRB_Stecker16_upper();
935 1 : em.setPhotonField(ebl);
936 1 : ebl = new IRB_Stecker16_lower();
937 1 : em.setPhotonField(ebl);
938 1 : ebl = new IRB_Finke22();
939 1 : em.setPhotonField(ebl);
940 1 : urb = new URB_Fixsen11();
941 1 : em.setPhotonField(urb);
942 1 : urb = new URB_Nitu21();
943 2 : em.setPhotonField(urb);
944 2 : }
945 :
946 1 : TEST(EMTripletPairProduction, limitNextStep) {
947 : // Test if the interaction limits the next propagation step.
948 1 : ref_ptr<PhotonField> cmb = new CMB();
949 1 : EMTripletPairProduction m(cmb);
950 1 : Candidate c(11, 1E17 * eV);
951 1 : c.setNextStep(std::numeric_limits<double>::max());
952 1 : m.process(&c);
953 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
954 2 : }
955 :
956 1 : TEST(EMTripletPairProduction, secondaries) {
957 : // Test if secondaries are correctly produced.
958 1 : ref_ptr<PhotonField> cmb = new CMB();
959 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
960 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
961 1 : EMPairProduction m(cmb);
962 1 : m.setHaveElectrons(true);
963 1 : m.setThinning(0.);
964 :
965 : std::vector<ref_ptr<PhotonField>> fields;
966 1 : fields.push_back(cmb);
967 1 : fields.push_back(irb);
968 1 : fields.push_back(urb);
969 :
970 : // loop over photon backgrounds
971 4 : for (int f = 0; f < fields.size(); f++) {
972 3 : m.setPhotonField(fields[f]);
973 :
974 : // loop over energies Ep = (1e9 - 1e23) eV
975 423 : for (int i = 0; i < 140; i++) {
976 :
977 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
978 420 : Candidate c(11, Ep);
979 420 : c.setCurrentStep(1e4 * Mpc); // use lower value so that the test can run faster
980 420 : m.process(&c);
981 :
982 : // pass if no interaction has occured (no tabulated rates)
983 420 : if (c.current.getEnergy() == Ep)
984 : continue;
985 :
986 : // expect positive energy of primary electron
987 0 : EXPECT_GT(c.current.getEnergy(), 0);
988 0 : double Etot = c.current.getEnergy();
989 :
990 : // expect electron / positron with energies 0 < E < Ephoton
991 0 : for (int j = 0; j < c.secondaries.size(); j++) {
992 0 : Candidate s = *c.secondaries[j];
993 0 : EXPECT_EQ(abs(s.current.getId()), 11);
994 0 : EXPECT_GT(s.current.getEnergy(), 0);
995 0 : EXPECT_LT(s.current.getEnergy(), Ep);
996 0 : Etot += s.current.getEnergy();
997 0 : }
998 :
999 : // test energy conservation
1000 0 : EXPECT_NEAR(Ep, Etot, 1e-9);
1001 420 : }
1002 : }
1003 2 : }
1004 :
1005 1 : TEST(EMTripletPairProduction, interactionTag) {
1006 1 : EMTripletPairProduction m(new CMB());
1007 :
1008 : // test default interactionTag
1009 1 : EXPECT_TRUE(m.getInteractionTag() == "EMTP");
1010 :
1011 : // test secondary tag
1012 1 : m.setHaveElectrons(true);
1013 2 : Candidate c(11, 1 * EeV);
1014 1 : m.performInteraction(&c);
1015 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMTP");
1016 :
1017 : // test custom tag
1018 1 : m.setInteractionTag("myTag");
1019 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
1020 1 : }
1021 :
1022 : // EMInverseComptonScattering -------------------------------------------------
1023 1 : TEST(EMInverseComptonScattering, allBackgrounds) {
1024 : // Test if interaction data files are loaded.
1025 1 : ref_ptr<PhotonField> cmb = new CMB();
1026 1 : EMInverseComptonScattering em(cmb);
1027 1 : ref_ptr<PhotonField> ebl = new IRB_Kneiske04();
1028 1 : em.setPhotonField(ebl);
1029 1 : ref_ptr<PhotonField> urb = new URB_Protheroe96();
1030 1 : em.setPhotonField(urb);
1031 1 : ebl = new IRB_Stecker05();
1032 1 : em.setPhotonField(ebl);
1033 1 : ebl = new IRB_Franceschini08();
1034 1 : em.setPhotonField(ebl);
1035 1 : ebl = new IRB_Finke10();
1036 1 : em.setPhotonField(ebl);
1037 1 : ebl = new IRB_Dominguez11();
1038 1 : em.setPhotonField(ebl);
1039 1 : ebl = new IRB_Gilmore12();
1040 1 : em.setPhotonField(ebl);
1041 1 : ebl = new IRB_Stecker16_upper();
1042 1 : em.setPhotonField(ebl);
1043 1 : ebl = new IRB_Stecker16_lower();
1044 1 : em.setPhotonField(ebl);
1045 1 : ebl = new IRB_Finke22();
1046 1 : em.setPhotonField(ebl);
1047 1 : urb = new URB_Fixsen11();
1048 1 : em.setPhotonField(urb);
1049 1 : urb = new URB_Nitu21();
1050 2 : em.setPhotonField(urb);
1051 2 : }
1052 :
1053 1 : TEST(EMInverseComptonScattering, limitNextStep) {
1054 : // Test if the interaction limits the next propagation step.
1055 1 : ref_ptr<PhotonField> cmb = new CMB();
1056 1 : EMInverseComptonScattering m(cmb);
1057 1 : Candidate c(11, 1E17 * eV);
1058 1 : c.setNextStep(std::numeric_limits<double>::max());
1059 1 : m.process(&c);
1060 1 : EXPECT_LT(c.getNextStep(), std::numeric_limits<double>::max());
1061 2 : }
1062 :
1063 1 : TEST(EMInverseComptonScattering, secondaries) {
1064 : // Test if secondaries are correctly produced.
1065 1 : ref_ptr<PhotonField> cmb = new CMB();
1066 1 : ref_ptr<PhotonField> irb = new IRB_Saldana21();
1067 1 : ref_ptr<PhotonField> urb = new URB_Nitu21();
1068 1 : EMPairProduction m(cmb);
1069 1 : m.setHaveElectrons(true);
1070 1 : m.setThinning(0.);
1071 :
1072 : std::vector<ref_ptr<PhotonField>> fields;
1073 1 : fields.push_back(cmb);
1074 1 : fields.push_back(irb);
1075 1 : fields.push_back(urb);
1076 :
1077 : // loop over photon backgrounds
1078 4 : for (int f = 0; f < fields.size(); f++) {
1079 3 : m.setPhotonField(fields[f]);
1080 :
1081 : // loop over energies Ep = (1e9 - 1e23) eV
1082 423 : for (int i = 0; i < 140; i++) {
1083 420 : double Ep = pow(10, 9.05 + 0.1 * i) * eV;
1084 420 : Candidate c(11, Ep);
1085 420 : c.setCurrentStep(1e3 * Mpc); // use lower value so that the test can run faster
1086 420 : m.process(&c);
1087 :
1088 : // pass if no interaction has occured (no tabulated rates)
1089 420 : if (c.current.getEnergy() == Ep)
1090 : continue;
1091 :
1092 : // expect positive energy of primary electron
1093 0 : EXPECT_GT(c.current.getEnergy(), 0);
1094 :
1095 : // expect photon with energy 0 < E < Ephoton
1096 0 : Candidate s = *c.secondaries[0];
1097 0 : EXPECT_EQ(abs(s.current.getId()), 22);
1098 0 : EXPECT_TRUE(s.current.getEnergy() >= 0.);
1099 0 : EXPECT_TRUE(s.current.getEnergy() < Ep);
1100 :
1101 :
1102 0 : double Etot = c.current.getEnergy();
1103 0 : for (int j = 0; j < c.secondaries.size(); j++) {
1104 0 : s = *c.secondaries[j];
1105 0 : Etot += s.current.getEnergy();
1106 : }
1107 0 : EXPECT_NEAR(Ep, Etot, 1e-9);
1108 420 : }
1109 : }
1110 2 : }
1111 :
1112 1 : TEST(EMInverseComptonScattering, interactionTag) {
1113 1 : EMInverseComptonScattering m(new CMB());
1114 :
1115 : // test default interactionTag
1116 1 : EXPECT_TRUE(m.getInteractionTag() == "EMIC");
1117 :
1118 : // test secondary tag
1119 1 : m.setHavePhotons(true);
1120 2 : Candidate c(11, 1 * PeV);
1121 1 : m.performInteraction(&c);
1122 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "EMIC");
1123 :
1124 : // test custom tag
1125 1 : m.setInteractionTag("myTag");
1126 1 : EXPECT_TRUE(m.getInteractionTag() == "myTag");
1127 1 : }
1128 :
1129 : // SynchrotronRadiation -------------------------------------------------
1130 1 : TEST(SynchrotronRadiation, interactionTag) {
1131 1 : SynchrotronRadiation s(1 * muG, true);
1132 :
1133 : // test default interactionTag
1134 1 : EXPECT_TRUE(s.getInteractionTag() == "SYN");
1135 :
1136 : // test secondary tag
1137 2 : Candidate c(11, 10 * PeV);
1138 1 : c.setCurrentStep(1 * pc);
1139 1 : s.process(&c);
1140 1 : EXPECT_TRUE(c.secondaries[0] -> getTagOrigin() == "SYN");
1141 :
1142 : // test custom tag
1143 1 : s.setInteractionTag("myTag");
1144 1 : EXPECT_TRUE(s.getInteractionTag() == "myTag");
1145 1 : }
1146 :
1147 1 : TEST(SynchrotronRadiation, simpleTestRMS) {
1148 : // test initialisation with B_rms
1149 :
1150 : // check default values
1151 1 : SynchrotronRadiation sync;
1152 :
1153 1 : EXPECT_EQ(sync.getBrms(), 0);
1154 1 : EXPECT_FALSE(sync.getHavePhotons());
1155 1 : EXPECT_EQ(sync.getThinning(), 0);
1156 1 : EXPECT_EQ(sync.getLimit(), 0.1);
1157 1 : EXPECT_EQ(sync.getMaximumSamples(), 0);
1158 1 : EXPECT_EQ(sync.getSecondaryThreshold(), 1 * MeV);
1159 :
1160 : // init with custom values
1161 1 : double b = 1 * muG;
1162 1 : double thinning = 0.23;
1163 1 : int samples = 4;
1164 1 : double limit = 0.123;
1165 2 : SynchrotronRadiation sync2(b, true, thinning, samples, limit);
1166 :
1167 1 : EXPECT_EQ(sync2.getBrms(), b);
1168 1 : EXPECT_TRUE(sync2.getHavePhotons());
1169 1 : EXPECT_EQ(sync2.getThinning(), thinning);
1170 1 : EXPECT_EQ(sync2.getLimit(), limit);
1171 1 : EXPECT_EQ(sync2.getMaximumSamples(), samples);
1172 1 : EXPECT_EQ(sync2.getSecondaryThreshold(), 1 * MeV);
1173 1 : }
1174 :
1175 1 : TEST(SynchrotronRadiation, simpleTestField) {
1176 : // test initialisation with field
1177 :
1178 : // check default values
1179 : Vector3d b(0, 0, 1 * muG);
1180 1 : ref_ptr<MagneticField> field = new UniformMagneticField(b);
1181 1 : SynchrotronRadiation sync(field);
1182 :
1183 1 : EXPECT_EQ(sync.getBrms(), 0);
1184 1 : EXPECT_FALSE(sync.getHavePhotons());
1185 1 : EXPECT_EQ(sync.getThinning(), 0);
1186 1 : EXPECT_EQ(sync.getLimit(), 0.1);
1187 1 : EXPECT_EQ(sync.getMaximumSamples(), 0);
1188 1 : EXPECT_EQ(sync.getSecondaryThreshold(), 1 * MeV);
1189 1 : Vector3d fieldAtPosition = sync.getField() -> getField(Vector3d(1, 2 , 3));
1190 1 : EXPECT_EQ(fieldAtPosition.getR(), b.getR());
1191 :
1192 : // init with custom values
1193 1 : double thinning = 0.23;
1194 1 : int samples = 4;
1195 1 : double limit = 0.123;
1196 2 : SynchrotronRadiation sync2(field, true, thinning, samples, limit);
1197 :
1198 1 : EXPECT_EQ(sync2.getBrms(), 0);
1199 1 : EXPECT_TRUE(sync2.getHavePhotons());
1200 1 : EXPECT_EQ(sync2.getThinning(), thinning);
1201 1 : EXPECT_EQ(sync2.getLimit(), limit);
1202 1 : EXPECT_EQ(sync2.getMaximumSamples(), samples);
1203 1 : EXPECT_EQ(sync2.getSecondaryThreshold(), 1 * MeV);
1204 1 : fieldAtPosition = sync2.getField() -> getField(Vector3d(1, 2 , 3));
1205 1 : EXPECT_EQ(fieldAtPosition.getR(), b.getR());
1206 2 : }
1207 :
1208 1 : TEST(SynchrotronRadiation, getSetFunctions) {
1209 1 : SynchrotronRadiation sync;
1210 :
1211 : // have photons
1212 1 : sync.setHavePhotons(true);
1213 1 : EXPECT_TRUE(sync.getHavePhotons());
1214 :
1215 : // Brms
1216 1 : sync.setBrms(5 * muG);
1217 1 : EXPECT_EQ(sync.getBrms(), 5 * muG);
1218 :
1219 : // thinning
1220 1 : sync.setThinning(0.345);
1221 1 : EXPECT_EQ(sync.getThinning(), 0.345);
1222 :
1223 : // limit
1224 1 : sync.setLimit(0.234);
1225 1 : EXPECT_EQ(sync.getLimit(), 0.234);
1226 :
1227 : // max samples
1228 1 : sync.setMaximumSamples(12345);
1229 1 : EXPECT_EQ(sync.getMaximumSamples(), 12345);
1230 :
1231 : // field
1232 : Vector3d b(1,2,3);
1233 1 : ref_ptr<MagneticField> field = new UniformMagneticField(b);
1234 1 : sync.setField(field);
1235 2 : EXPECT_TRUE(field == sync.getField()); // same pointer
1236 :
1237 : // secondary threshold
1238 1 : sync.setSecondaryThreshold(1 * eV);
1239 1 : EXPECT_EQ(sync.getSecondaryThreshold(), 1 * eV);
1240 1 : }
1241 :
1242 1 : TEST(SynchrotronRadiation, energyLoss) {
1243 : double brms = 1 * muG;
1244 : double step = 1 * kpc;
1245 1 : SynchrotronRadiation sync(brms, false);
1246 :
1247 : double dE, lf, Rg, dEdx;
1248 1 : Candidate c(11);
1249 1 : c.setCurrentStep(step);
1250 1 : c.setNextStep(step);
1251 : double charge = eplus;
1252 :
1253 : // 1 GeV
1254 1 : c.current.setEnergy(1 * GeV);
1255 1 : lf = c.current.getLorentzFactor();
1256 : Rg = 1 * GeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1257 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1258 1 : dE = dEdx * step;
1259 1 : sync.process(&c);
1260 1 : EXPECT_NEAR(1 * GeV - c.current.getEnergy(), dE, 0.01 * dE);
1261 :
1262 : // 100 GeV
1263 1 : c.current.setEnergy(100 * GeV);
1264 1 : lf = c.current.getLorentzFactor();
1265 : Rg = 100 * GeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1266 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1267 1 : dE = dEdx * step;
1268 1 : sync.process(&c);
1269 1 : EXPECT_NEAR(100 * GeV - c.current.getEnergy(), dE, 0.01 * dE);
1270 :
1271 : // 10 TeV
1272 1 : c.current.setEnergy(10 * TeV);
1273 1 : lf = c.current.getLorentzFactor();
1274 : Rg = 10 * TeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1275 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1276 1 : dE = dEdx * step;
1277 1 : sync.process(&c);
1278 1 : EXPECT_NEAR(10 * TeV - c.current.getEnergy(), dE, 0.01 * dE);
1279 :
1280 : // 1 PeV
1281 1 : c.current.setEnergy(1 * PeV);
1282 1 : lf = c.current.getLorentzFactor();
1283 : Rg = 1 * PeV / charge / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1284 1 : dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
1285 1 : dE = dEdx * step;
1286 1 : sync.process(&c);
1287 1 : EXPECT_NEAR(1 * PeV - c.current.getEnergy(), dE, 0.01 * dE);
1288 1 : }
1289 :
1290 1 : TEST(SynchrotronRadiation, PhotonEnergy) {
1291 : double brms = 1 * muG;
1292 1 : SynchrotronRadiation sync(brms, true);
1293 1 : sync.setSecondaryThreshold(0.); // allow all secondaries for testing
1294 1 : sync.setMaximumSamples(1000); // reduce the amount of generated secondaries
1295 :
1296 : double E = 1 * TeV;
1297 1 : Candidate c(11, E);
1298 1 : c.setCurrentStep(10 * pc);
1299 1 : c.setNextStep(10 * pc);
1300 :
1301 1 : double lf = c.current.getLorentzFactor();
1302 : double Rg = E / eplus / c_light / (brms * sqrt(2. / 3) ); // factor 2/3 for avg magnetic field direction.
1303 1 : double Ecrit = 3. / 4 * h_planck / M_PI * c_light * pow(lf, 3) / Rg;
1304 :
1305 1 : sync.process(&c);
1306 1 : EXPECT_TRUE(c.secondaries.size() > 0); // must have secondaries
1307 :
1308 : // check avg energy of the secondary photons
1309 : double Esec = 0;
1310 : double weightSum = 0;
1311 1001 : for (size_t i = 0; i < c.secondaries.size(); i++) {
1312 1000 : double weight = c.secondaries[i]->getWeight();
1313 1000 : Esec += c.secondaries[i] -> current.getEnergy()*weight;
1314 1000 : weightSum += weight;
1315 : }
1316 1 : Esec /= weightSum;
1317 :
1318 1 : EXPECT_NEAR(Esec, Ecrit, Ecrit);
1319 1 : }
1320 :
1321 0 : int main(int argc, char **argv) {
1322 0 : ::testing::InitGoogleTest(&argc, argv);
1323 0 : return RUN_ALL_TESTS();
1324 : }
1325 :
1326 : } // namespace crpropa
|
