Line data Source code
1 : /// Tests gamale with a test lense from file and demonstrating basic
2 : /// usage
3 :
4 : //--------------------------------------------
5 : // Project: Galactic magnetic Lense (GaMaLe) -
6 : // Copyright (C) 2009 Tobias Winchen -
7 : // RWTH Aachen, Germany -
8 : // Contact: winchen@physik.rwth-achen.de -
9 : // Licensed under the GNU GPL v2 -
10 : //--------------------------------------------
11 :
12 : #include <string>
13 : #include "gtest/gtest.h"
14 :
15 : #include "crpropa/magneticLens/MagneticLens.h"
16 : #include "crpropa/magneticLens/ModelMatrix.h"
17 : #include "crpropa/magneticLens/Pixelization.h"
18 : #include "crpropa/magneticLens/ParticleMapsContainer.h"
19 : #include "crpropa/Common.h"
20 :
21 : using namespace std;
22 : using namespace crpropa;
23 :
24 1 : TEST(MagneticLens, Deflection)
25 : {
26 1 : MagneticLens magneticLens(5);
27 1 : Pixelization P(5);
28 1 : ModelMatrixType M;
29 1 : M.resize(P.nPix(), P.nPix());
30 1 : M.reserve(P.nPix());
31 :
32 : // Map any direction (p,t) to (p, -t)
33 12289 : for (int i=0;i<P.nPix();i++)
34 : {
35 : double theta, phi;
36 12288 : P.pix2Direction(i, phi, theta);
37 12288 : theta*= -1;
38 12288 : int j = P.direction2Pix(phi, theta);
39 12288 : M.insert(i,j) =1;
40 : }
41 :
42 1 : magneticLens.setLensPart(M, 10 * EeV, 100 * EeV);
43 :
44 12289 : for (int i=0; i < P.nPix(); i++)
45 : {
46 : double theta, phi;
47 12288 : P.pix2Direction(i, phi, theta);
48 12288 : double theta0 = theta;
49 :
50 : // No CR is allowed to be lost in this lens
51 12288 : EXPECT_TRUE(magneticLens.transformCosmicRay(20 * EeV, phi, theta));
52 12288 : EXPECT_NEAR(theta+theta0,0., 0.05);
53 : }
54 2 : }
55 :
56 1 : TEST(MagneticLens, Vector3Deflection)
57 : {
58 1 : MagneticLens magneticLens(5);
59 1 : Pixelization P(5);
60 1 : ModelMatrixType M;
61 1 : M.resize(P.nPix(), P.nPix());
62 1 : M.reserve(P.nPix());
63 :
64 : // No deflection
65 12289 : for (int i=0;i<P.nPix();i++)
66 : {
67 12288 : M.insert(i,i) = 1;
68 : }
69 :
70 1 : magneticLens.setLensPart(M, 10 * EeV, 100 * EeV);
71 :
72 : Vector3d u(1,0,0);
73 : Vector3d v(u);
74 1 : EXPECT_TRUE(magneticLens.transformCosmicRay(20 * EeV, v));
75 1 : EXPECT_NEAR(u.getAngleTo(v), 0., 2. / 180 * M_PI);
76 :
77 1 : u.x = 0; u.y = 1; u.z = 0; v = u;
78 1 : EXPECT_TRUE(magneticLens.transformCosmicRay(20 * EeV, v));
79 1 : EXPECT_NEAR(u.getAngleTo(v), 0., 2. / 180 * M_PI);
80 :
81 1 : u.x = 0; u.y = 0; u.z = 1; v = u;
82 1 : EXPECT_TRUE(magneticLens.transformCosmicRay(20 * EeV, v));
83 1 : EXPECT_NEAR(u.getAngleTo(v), 0., 2. / 180 * M_PI);
84 :
85 1 : u.x = 1; u.y = -2; u.z = 3; v = u;
86 1 : EXPECT_TRUE(magneticLens.transformCosmicRay(20 * EeV, v));
87 1 : EXPECT_NEAR(u.getAngleTo(v), 0., 2. / 180 * M_PI);
88 2 : }
89 :
90 1 : TEST(MagneticLens, OutOfBoundsEnergy)
91 : {
92 1 : MagneticLens magneticLens(5);
93 1 : Pixelization P(5);
94 1 : ModelMatrixType M;
95 1 : M.resize(P.nPix(), P.nPix());
96 1 : M.reserve(P.nPix());
97 1 : magneticLens.setLensPart(M,10. * EeV, 100. * EeV);
98 1 : double theta = 0, phi = 0;
99 1 : EXPECT_FALSE(magneticLens.transformCosmicRay(1. * EeV, phi, theta));
100 2 : }
101 :
102 :
103 1 : TEST(Pixelization, angularDistance)
104 : {
105 : // test for correct angular distance in case of same vectors
106 1 : Pixelization P(6);
107 49153 : for (int idx =0; idx < P.nPix(); idx++)
108 : {
109 49152 : double ang = P.angularDistance(idx,idx);
110 49152 : EXPECT_TRUE(ang == ang);
111 : }
112 1 : }
113 :
114 :
115 1 : TEST(ParticleMapsContainer, addParticle)
116 : {
117 1 : ParticleMapsContainer maps;
118 1 : maps.addParticle(1000010010, 1 * EeV, 0 , 0 );
119 1 : std::vector<int> pids = maps.getParticleIds();
120 1 : EXPECT_EQ(pids.size(), 1);
121 1 : EXPECT_EQ(pids[0], 1000010010);
122 :
123 1 : std::vector<double> energies = maps.getEnergies(1000010010);
124 1 : EXPECT_EQ(energies.size(), 1);
125 1 : }
126 :
127 1 : TEST(ParticleMapsContainer, getRandomParticles)
128 : {
129 1 : ParticleMapsContainer maps(0.002);
130 1 : maps.addParticle(1000010010, 1 * EeV, 0 , 0 );
131 : std::vector<double> energies;
132 : std::vector<double> lons;
133 : std::vector<double> lats;
134 : std::vector<int> particleIds;
135 :
136 1 : size_t N = 420;
137 1 : maps.getRandomParticles(N, particleIds, energies, lons, lats);
138 :
139 1 : EXPECT_EQ(energies.size(), N);
140 1 : EXPECT_EQ(lons.size(), N);
141 1 : EXPECT_EQ(lats.size(), N);
142 1 : EXPECT_EQ(particleIds.size(), N);
143 :
144 421 : for(size_t i = 0; i < N; i++)
145 : {
146 420 : EXPECT_NEAR(log10(energies[i]), 18, 0.002);
147 420 : EXPECT_EQ(particleIds[i], 1000010010);
148 420 : EXPECT_NEAR(lons[i], 0, 2./180*M_PI);
149 420 : EXPECT_NEAR(lats[i], 0, 2./180*M_PI);
150 : }
151 :
152 1 : }
153 :
154 1 : TEST(Pixelization, randomDirectionInPixel)
155 : {
156 1 : Pixelization p(6);
157 : const double long0 = -35./180 * M_PI;
158 : const double lat0 = 12.08/180 * M_PI;
159 :
160 1 : int pix = p.direction2Pix(long0, lat0);
161 :
162 : double rlon, rlat;
163 1 : p.getRandomDirectionInPixel(pix, rlon, rlat);
164 :
165 : // new direction should be inside the pixel
166 1 : EXPECT_EQ(pix, p.direction2Pix(rlon, rlat));
167 :
168 : // new direction should be different from input
169 1 : EXPECT_FALSE(long0 == rlon);
170 1 : EXPECT_FALSE(lat0 == rlat);
171 :
172 1 : }
|
