Line data Source code
1 : #include "crpropa/Source.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/ParticleID.h"
4 :
5 : #include "gtest/gtest.h"
6 : #include <stdexcept>
7 :
8 : namespace crpropa {
9 :
10 2 : TEST(SourcePosition, simpleTest) {
11 : Vector3d position(1, 2, 3);
12 1 : SourcePosition source(position);
13 1 : ParticleState ps;
14 1 : source.prepareParticle(ps);
15 2 : EXPECT_EQ(position, ps.getPosition());
16 1 : }
17 :
18 1 : TEST(SourceMultiplePositions, simpleTest) {
19 1 : SourceMultiplePositions source;
20 1 : source.add(Vector3d(1, 0, 0), 0.25);
21 1 : source.add(Vector3d(2, 0, 0), 0.75);
22 1 : ParticleState ps;
23 : int n1 = 0;
24 : int n2 = 0;
25 10001 : for (int i = 0; i < 10000; i++) {
26 10000 : source.prepareParticle(ps);
27 10000 : if (ps.getPosition().x == 1)
28 2536 : n1++;
29 7464 : else if (ps.getPosition().x == 2)
30 7464 : n2++;
31 : }
32 1 : EXPECT_NEAR(n1, 2500, 5 * sqrt(2500));
33 1 : EXPECT_NEAR(n2, 7500, 5 * sqrt(7500));
34 1 : }
35 :
36 2 : TEST(SourceUniformSphere, simpleTest) {
37 : Vector3d center(0, 0, 0);
38 1 : double radius = 110;
39 1 : SourceUniformSphere source(center, radius);
40 1 : ParticleState ps;
41 1 : source.prepareParticle(ps);
42 1 : double distance = ps.getPosition().getDistanceTo(center);
43 1 : EXPECT_GE(radius, distance);
44 1 : }
45 :
46 2 : TEST(SourceUniformHollowSphere, simpleTest) {
47 : Vector3d center(0, 0, 0);
48 1 : double radius_inner = 50;
49 1 : double radius_outer = 110;
50 : SourceUniformHollowSphere source(center,
51 : radius_inner,
52 1 : radius_outer);
53 101 : for (int i=0; i < 100; ++i) {
54 100 : ParticleState ps;
55 100 : source.prepareParticle(ps);
56 100 : double distance = ps.getPosition().getDistanceTo(center);
57 100 : EXPECT_GE(radius_outer, distance);
58 100 : EXPECT_LE(radius_inner, distance);
59 : }
60 1 : }
61 :
62 2 : TEST(SourceUniformBox, simpleTest) {
63 : Vector3d origin(-7, -2, 0);
64 : Vector3d size(13, 55, 192);
65 1 : SourceUniformBox box(origin, size);
66 1 : ParticleState p;
67 1 : box.prepareParticle(p);
68 1 : Vector3d pos = p.getPosition();
69 1 : EXPECT_LE(origin.x, pos.x);
70 1 : EXPECT_LE(origin.y, pos.y);
71 1 : EXPECT_LE(origin.z, pos.z);
72 1 : EXPECT_GE(size.x, pos.x);
73 1 : EXPECT_GE(size.y, pos.y);
74 1 : EXPECT_GE(size.z, pos.z);
75 1 : }
76 :
77 2 : TEST(SourceUniformCylinder, simpleTest) {
78 : Vector3d center(0, 0, 0);
79 : double radius = 15;
80 : double height = 2;
81 1 : SourceUniformCylinder cylinder(center, height, radius);
82 1 : ParticleState ps;
83 1 : cylinder.prepareParticle(ps);
84 1 : Vector3d pos = ps.getPosition();
85 1 : double R2 = pos.x*pos.x+pos.y*pos.y;
86 1 : double H = pow(pos.z*pos.z, 0.5);
87 1 : EXPECT_GE(radius*radius, R2);
88 1 : EXPECT_GE(height/2., H);
89 1 : }
90 :
91 1 : TEST(SourceSNRDistribution, simpleTest) {
92 : double R_earth = 8.5*kpc;
93 : double alpha = 2.0;
94 : double beta = 3.53;
95 : double Z_G = 0.3*kpc;
96 1 : SourceSNRDistribution snr(R_earth,alpha, beta, Z_G);
97 1 : ParticleState ps;
98 1 : snr.prepareParticle(ps);
99 1 : Vector3d pos = ps.getPosition();
100 1 : double R2 = pos.x*pos.x+pos.y*pos.y;
101 1 : EXPECT_GE(20*kpc*20*kpc, R2); // radius must be smaller than 20 kpc
102 :
103 : double R2_mean = 0.;
104 : double Z_mean = 0.;
105 100001 : for (size_t i=0; i<100000; i++) {
106 100000 : snr.prepareParticle(ps);
107 100000 : Vector3d pos = ps.getPosition();
108 100000 : R2_mean += pow(pos.x/kpc, 2.)+pow(pos.y/kpc, 2.);
109 100000 : Z_mean += pos.z/kpc;
110 : }
111 1 : R2_mean/=100000.;
112 1 : Z_mean/=100000.;
113 1 : EXPECT_NEAR(100.306, R2_mean, 1);
114 1 : EXPECT_NEAR(0., Z_mean, 0.1);
115 1 : }
116 :
117 2 : TEST(SourceDensityGrid, withInRange) {
118 : // Create a grid with 10^3 cells ranging from (0, 0, 0) to (10, 10, 10)
119 : Vector3d origin(0, 0, 0);
120 : int cells = 10;
121 : double spacing = 1;
122 1 : auto grid = new Grid1f(origin, cells, spacing);
123 11 : for (int ix = 0; ix < cells; ix++)
124 110 : for (int iy = 0; iy < cells; iy++)
125 1100 : for (int iz = 0; iz < cells; iz++)
126 1000 : grid->get(ix, iy, iz) = ix * iy * iz;
127 :
128 2 : SourceDensityGrid source(grid);
129 1 : ParticleState p;
130 :
131 1 : source.prepareParticle(p);
132 1 : Vector3d pos = p.getPosition();
133 :
134 : // dialed positions should be within the volume (0, 0, 0) - (10, 10, 10)
135 1 : EXPECT_LE(0, pos.x);
136 1 : EXPECT_GE(10, pos.x);
137 1 : EXPECT_LE(0, pos.y);
138 1 : EXPECT_GE(10, pos.y);
139 1 : EXPECT_LE(0, pos.z);
140 1 : EXPECT_GE(10, pos.z);
141 1 : }
142 :
143 2 : TEST(SourceDensityGrid, OneAllowedCell) {
144 : // Create a grid with 2^3 cells ranging from (0, 0, 0) to (4, 4, 4)
145 : Vector3d origin(0, 0, 0);
146 : int cells = 2;
147 : double spacing = 2;
148 1 : auto grid = new Grid1f(origin, cells, spacing);
149 :
150 : // set all but one cells to 0
151 3 : for (int ix = 0; ix < cells; ix++)
152 6 : for (int iy = 0; iy < cells; iy++)
153 12 : for (int iz = 0; iz < cells; iz++)
154 8 : grid->get(ix, iy, iz) = 0;
155 :
156 : // set the first cell ((0, 0, 0) to (2, 2, 2))
157 1 : grid->get(0, 0, 0) = 1;
158 :
159 2 : SourceDensityGrid source(grid);
160 1 : ParticleState p;
161 :
162 1 : int nFalse = 0;
163 : Vector3d mean(0, 0, 0);
164 10001 : for (int i = 0; i < 10000; i++) {
165 10000 : source.prepareParticle(p);
166 10000 : Vector3d pos = p.getPosition();
167 : mean += pos;
168 10000 : if ((pos.x < 0) or (pos.x > 2) or (pos.y < 0) or (pos.y > 2)
169 10000 : or (pos.z < 0) or (pos.z > 2))
170 0 : nFalse++;
171 : }
172 :
173 : // only the first bin should get dialed
174 1 : EXPECT_EQ(0, nFalse);
175 :
176 : // mean should be close to (1, 1, 1) if random positions are uniform in (0, 0, 0) - (2, 2, 2)
177 : mean /= 10000;
178 1 : EXPECT_NEAR(1, mean.x, 0.2);
179 1 : EXPECT_NEAR(1, mean.y, 0.2);
180 1 : EXPECT_NEAR(1, mean.z, 0.2);
181 1 : }
182 :
183 2 : TEST(SourceDensityGrid1D, withInRange) {
184 : // Create a grid with 10 cells ranging from 0 to 10
185 : Vector3d origin(0, 0, 0);
186 : int nCells = 10;
187 : double spacing = 1.;
188 1 : auto grid = new Grid1f(origin, nCells, 1, 1, spacing);
189 :
190 : // set some values
191 11 : for (int i = 0; i < 10; i++) {
192 10 : grid->get(i, 0, 0) = 2;
193 : }
194 :
195 2 : SourceDensityGrid1D source(grid);
196 1 : ParticleState p;
197 :
198 1 : source.prepareParticle(p);
199 1 : Vector3d pos = p.getPosition();
200 : // dialed position should be within the range 0 - 10
201 1 : EXPECT_LE(0, pos.x);
202 1 : EXPECT_GE(10, pos.x);
203 1 : }
204 :
205 2 : TEST(SourceDensityGrid1D, OneAllowedCell) {
206 : // Test if the only allowed cells is repeatedly selected
207 : Vector3d origin(0, 0, 0);
208 : int nCells = 10;
209 : double spacing = 1.;
210 1 : auto grid = new Grid1f(origin, nCells, 1, 1, spacing);
211 :
212 : // set some values
213 11 : for (int i = 0; i < 10; i++) {
214 10 : grid->get(i, 0, 0) = 0;
215 : }
216 1 : grid->get(5, 0, 0) = 1;
217 :
218 2 : SourceDensityGrid1D source(grid);
219 1 : ParticleState p;
220 :
221 101 : for (int i = 0; i < 100; i++) {
222 100 : source.prepareParticle(p);
223 : // dialed position should be in range 5-6
224 100 : Vector3d pos = p.getPosition();
225 100 : EXPECT_LE(5, pos.x);
226 100 : EXPECT_GE(6, pos.x);
227 : }
228 1 : }
229 :
230 1 : TEST(SourcePowerLawSpectrum, simpleTest) {
231 1 : double Emin = 4 * EeV;
232 1 : double Emax = 200 * EeV;
233 : double index = -2.7;
234 1 : SourcePowerLawSpectrum spectrum(Emin, Emax, index);
235 1 : ParticleState ps;
236 1 : spectrum.prepareParticle(ps);
237 :
238 : // energy should be within Emin - Emax
239 1 : EXPECT_LE(Emin, ps.getEnergy());
240 1 : EXPECT_GE(Emax, ps.getEnergy());
241 1 : }
242 :
243 1 : TEST(SourceComposition, simpleTest) {
244 1 : double Emin = 10;
245 : double Rmax = 100;
246 : double index = -1;
247 1 : SourceComposition source(Emin, Rmax, index);
248 1 : source.add(nucleusId(6, 3), 1);
249 1 : ParticleState p;
250 1 : source.prepareParticle(p);
251 1 : EXPECT_EQ(nucleusId(6, 3), p.getId());
252 1 : EXPECT_LE(Emin, p.getEnergy());
253 1 : EXPECT_GE(6 * Rmax, p.getEnergy());
254 1 : }
255 :
256 2 : TEST(SourceDirectedEmission, simpleTest) {
257 : Vector3d mu(1., 0., 0.);
258 : double kappa = 1000.;
259 1 : SourceDirectedEmission source(mu, kappa);
260 1 : Candidate c;
261 : Vector3d meanDir(0., 0., 0.);
262 1001 : for (size_t i = 0; i < 1000; i++) {
263 1000 : source.prepareCandidate(c);
264 1000 : meanDir += c.source.getDirection();
265 1000 : double w = c.getWeight();
266 1000 : EXPECT_GE(w, 0.);
267 : }
268 : meanDir /= 1000.;
269 1 : EXPECT_NEAR(meanDir.x, 1., 0.1);
270 1 : EXPECT_NEAR(meanDir.y, 0., 0.01);
271 1 : EXPECT_NEAR(meanDir.z, 0., 0.01);
272 2 : }
273 :
274 2 : TEST(SourceEmissionCone, simpleTest) {
275 : Vector3d direction(42., 0., 0.);
276 1 : double aperture = 1/42.;
277 :
278 1 : SourceEmissionCone source(direction, aperture);
279 :
280 1 : ParticleState p;
281 1 : source.prepareParticle(p);
282 1 : double angle = direction.getAngleTo(p.getDirection());
283 1 : EXPECT_LE(angle, aperture);
284 1 : }
285 :
286 : #ifdef CRPROPA_HAVE_MUPARSER
287 1 : TEST(SourceGenericComposition, simpleTest) {
288 : double Emin = 10;
289 : double Emax = 100;
290 1 : SourceGenericComposition source(Emin, Emax, "E^-2");
291 1 : int id1 = nucleusId(6, 3);
292 1 : int id2 = nucleusId(12, 6);
293 1 : source.add(id1, 1);
294 1 : source.add(id2, 10);
295 1 : ParticleState p;
296 : size_t id1Count = 0, id2Count = 0;
297 : size_t ElowCount = 0, EhighCount = 0;
298 1 : size_t n = 100000;
299 100001 : for (size_t i = 0; i < n; i++) {
300 100000 : source.prepareParticle(p);
301 100000 : if (p.getId() == id1)
302 8937 : id1Count++;
303 100000 : if (p.getId() == id2)
304 91063 : id2Count++;
305 100000 : double e = p.getEnergy();
306 100000 : if ( (e >= Emin) && (e < 20))
307 55771 : ElowCount++;
308 100000 : if ( (e >= 20) && (e <= Emax))
309 44229 : EhighCount++;
310 :
311 : }
312 1 : EXPECT_EQ(n, id1Count + id2Count);
313 1 : EXPECT_EQ(n, ElowCount + EhighCount);
314 1 : EXPECT_NEAR((float)id1Count/(float)id2Count, 0.1, 0.01);
315 1 : EXPECT_NEAR((float)ElowCount/(float)EhighCount, 1.25, 0.1);
316 1 : }
317 : #endif
318 :
319 1 : TEST(SourceComposition, throwNoIsotope) {
320 1 : SourceComposition source(1, 10, -1);
321 1 : ParticleState ps;
322 1 : EXPECT_THROW(source.prepareParticle(ps), std::runtime_error);
323 1 : }
324 :
325 1 : TEST(SourceRedshiftEvolution, testInRange) {
326 1 : Candidate c;
327 :
328 1 : double zmin = 0.5;
329 1 : double zmax = 2.5;
330 :
331 : // general case: m
332 1 : SourceRedshiftEvolution source1(3.2, zmin, zmax);
333 101 : for (int i = 0; i < 100; i++) {
334 100 : source1.prepareCandidate(c);
335 100 : EXPECT_LE(zmin, c.getRedshift());
336 100 : EXPECT_GE(zmax, c.getRedshift());
337 : }
338 :
339 : // general case: m = -1
340 1 : SourceRedshiftEvolution source2(-1, zmin, zmax);
341 101 : for (int i = 0; i < 100; i++) {
342 100 : source2.prepareCandidate(c);
343 100 : EXPECT_LE(zmin, c.getRedshift());
344 100 : EXPECT_GE(zmax, c.getRedshift());
345 : }
346 1 : }
347 :
348 2 : TEST(Source, allPropertiesUsed) {
349 : Source source;
350 1 : source.add(new SourcePosition(Vector3d(10, 0, 0) * Mpc));
351 2 : source.add(new SourceIsotropicEmission());
352 1 : source.add(new SourcePowerLawSpectrum(5 * EeV, 100 * EeV, -2));
353 1 : source.add(new SourceParticleType(nucleusId(8, 4)));
354 1 : source.add(new SourceRedshift(2));
355 :
356 1 : Candidate c = *source.getCandidate();
357 :
358 1 : EXPECT_EQ(2, c.getRedshift());
359 :
360 1 : ParticleState p;
361 :
362 : p = c.source;
363 1 : EXPECT_EQ(nucleusId(8, 4), p.getId());
364 1 : EXPECT_LE(5 * EeV, p.getEnergy());
365 1 : EXPECT_GE(100 * EeV, p.getEnergy());
366 2 : EXPECT_EQ(Vector3d(10, 0, 0) * Mpc, p.getPosition());
367 :
368 : p = c.created;
369 1 : EXPECT_EQ(nucleusId(8, 4), p.getId());
370 1 : EXPECT_LE(5 * EeV, p.getEnergy());
371 1 : EXPECT_GE(100 * EeV, p.getEnergy());
372 2 : EXPECT_EQ(Vector3d(10, 0, 0) * Mpc, p.getPosition());
373 :
374 : p = c.previous;
375 1 : EXPECT_EQ(nucleusId(8, 4), p.getId());
376 1 : EXPECT_LE(5 * EeV, p.getEnergy());
377 1 : EXPECT_GE(100 * EeV, p.getEnergy());
378 2 : EXPECT_EQ(Vector3d(10, 0, 0) * Mpc, p.getPosition());
379 :
380 : p = c.current;
381 1 : EXPECT_EQ(nucleusId(8, 4), p.getId());
382 1 : EXPECT_LE(5 * EeV, p.getEnergy());
383 1 : EXPECT_GE(100 * EeV, p.getEnergy());
384 2 : EXPECT_EQ(Vector3d(10, 0, 0) * Mpc, p.getPosition());
385 2 : }
386 :
387 2 : TEST(SourceList, simpleTest) {
388 : // test if source list works with one source
389 : SourceList sourceList;
390 1 : ref_ptr<Source> source = new Source;
391 1 : source->add(new SourcePosition(Vector3d(10, 0, 0)));
392 1 : sourceList.add(source);
393 :
394 1 : ref_ptr<Candidate> c = sourceList.getCandidate();
395 :
396 2 : EXPECT_EQ(Vector3d(10, 0, 0), c->created.getPosition());
397 2 : EXPECT_EQ(Vector3d(10, 0, 0), c->previous.getPosition());
398 2 : EXPECT_EQ(Vector3d(10, 0, 0), c->current.getPosition());
399 1 : }
400 :
401 2 : TEST(SourceList, noSource) {
402 : // test if an error is thrown when source list empty
403 : SourceList sourceList;
404 1 : EXPECT_THROW(sourceList.getCandidate(), std::runtime_error);
405 1 : }
406 :
407 2 : TEST(SourceList, luminosity) {
408 : // test if the sources are dialed according to their luminosities
409 : SourceList sourceList;
410 :
411 1 : ref_ptr<Source> source1 = new Source;
412 1 : source1->add(new SourceEnergy(100));
413 1 : sourceList.add(source1, 80);
414 :
415 1 : ref_ptr<Source> source2 = new Source;
416 1 : source2->add(new SourceEnergy(0));
417 1 : sourceList.add(source2, 20);
418 :
419 : double meanE = 0;
420 1001 : for (int i = 0; i < 1000; i++) {
421 1000 : ref_ptr<Candidate> c = sourceList.getCandidate();
422 1000 : meanE += c->created.getEnergy();
423 : }
424 1 : meanE /= 1000;
425 1 : EXPECT_NEAR(80, meanE, 4); // this test can stochastically fail
426 1 : }
427 :
428 1 : TEST(SourceTag, sourceTag) {
429 1 : SourceTag tag("mySourceTag");
430 1 : Candidate c;
431 1 : tag.prepareCandidate(c);
432 2 : EXPECT_TRUE(c.getTagOrigin() == "mySourceTag");
433 1 : }
434 :
435 0 : int main(int argc, char **argv) {
436 0 : ::testing::InitGoogleTest(&argc, argv);
437 0 : return RUN_ALL_TESTS();
438 : }
439 :
440 : } // namespace crpropa
|
