Line data Source code
1 : //----------------------------------------------------------------------
2 : // This file is part of PARSEC (http://physik.rwth-aachen.de/parsec)
3 : // a parametrized simulation engine for cosmic rays.
4 : //
5 : // Copyright (C) 2011 Martin Erdmann, Peter Schiffer, Tobias Winchen
6 : // RWTH Aachen University, Germany
7 : // Contact: winchen@physik.rwth-aachen.de
8 : //
9 : // This program is free software: you can redistribute it and/or
10 : // modify it under the terms of the GNU General Public License as
11 : // published by the Free Software Foundation, either version 3 of
12 : // the License, or (at your option) any later version.
13 : //
14 : // This program is distributed in the hope that it will be useful,
15 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : // GNU General Public License for more details.
18 : //
19 : // You should have received a copy of the GNU General Public License
20 : // along with this program. If not, see <http://www.gnu.org/licenses/>.
21 : //----------------------------------------------------------------------
22 :
23 : #ifndef MAGNETICLENS_HH
24 : #define MAGNETICLENS_HH
25 :
26 : #include "crpropa/magneticLens/ModelMatrix.h"
27 : #include "crpropa/magneticLens/Pixelization.h"
28 : #include "crpropa/Units.h"
29 : #include "crpropa/Vector3.h"
30 :
31 : #include <vector>
32 : #include <string>
33 : #include <stdexcept>
34 : #include <fstream>
35 : #include <sstream>
36 : #include <iostream>
37 : #include <stdint.h>
38 :
39 :
40 : namespace crpropa
41 : {
42 :
43 : /// Holds one matrix for the lens and information about the rigidity range
44 : class LensPart
45 : {
46 : string _filename;
47 : double _rigidityMin;
48 : double _rigidityMax;
49 : ModelMatrixType M;
50 : double _maximumSumOfColumns;
51 : bool _maximumSumOfColumns_calculated;
52 :
53 : public:
54 0 : LensPart()
55 0 : {
56 0 : }
57 : /// File containing the matrix to be used in the range rigidityMin,
58 : /// rigidityMax in Joule
59 3 : LensPart(const std::string &filename, double rigidityMin, double rigidityMax) :
60 3 : _filename(filename), _rigidityMin(rigidityMin), _rigidityMax(rigidityMax), _maximumSumOfColumns_calculated(
61 3 : false), _maximumSumOfColumns(0)
62 : {
63 3 : }
64 :
65 3 : ~LensPart()
66 : {
67 3 : }
68 :
69 : /// Loads the matrix from file
70 : void loadMatrixFromFile()
71 : {
72 0 : deserialize(_filename, M);
73 0 : }
74 :
75 : /// Returns the filename of the matrix
76 : const std::string& getFilename()
77 : {
78 : return _filename;
79 : }
80 :
81 : /// Calculates the maximum of the sums of columns for the matrix
82 : double getMaximumOfSumsOfColumns()
83 : {
84 0 : if (!_maximumSumOfColumns_calculated)
85 : { // lazy calculation of maximum
86 0 : _maximumSumOfColumns = maximumOfSumsOfColumns(M);
87 0 : _maximumSumOfColumns_calculated = true;
88 : }
89 0 : return _maximumSumOfColumns;
90 : }
91 :
92 : /// Returns the minimum of the rigidity range for the lenspart in eV
93 : double getMinimumRigidity()
94 : {
95 12293 : return _rigidityMin / eV;
96 : }
97 :
98 : /// Returns the maximum of the rigidity range for the lenspart in eV
99 : double getMaximumRigidity()
100 : {
101 12292 : return _rigidityMax / eV;
102 : }
103 :
104 : /// Returns the modelmatrix
105 : ModelMatrixType& getMatrix()
106 : {
107 0 : return M;
108 : }
109 :
110 : /// Sets the modelmatrix
111 : void setMatrix(const ModelMatrixType& m)
112 : {
113 3 : M = m;
114 0 : }
115 :
116 :
117 : };
118 :
119 : /// Function to calculate the mean deflection [rad] of the matrix M, given a pixelization
120 : //double calculateMeanDeflection(const ModelMatrix &M,
121 : // const Pixelization &pixelization)
122 : //{
123 : // double totalDeflection = 0;
124 : // double weightSum = 0;
125 : // for (const_i2_t it1 = M.begin2(); it1 != (M.end2()); it1++)
126 : // {
127 : // for (const_i1_t it2 = it1.begin(); it2 != it1.end(); it2++)
128 : // {
129 : // totalDeflection+= pixelization.angularDistance(it2.index1(),
130 : // it2.index2()) * (*it2) ;
131 : // weightSum+= (*it2);
132 : // }
133 : // }
134 : // return totalDeflection / weightSum;
135 : //}
136 :
137 : typedef std::vector<LensPart*>::iterator LensPartIter;
138 : typedef std::vector<LensPart*>::const_iterator const_LensPartIter;
139 :
140 : /**
141 : * \addtogroup MagneticLenses
142 : * @{
143 : */
144 :
145 : /// The lens for the galactic magnetic field.
146 : /// Note that the energies refer to protons (Z=1). To be used with other particles with a different charge number please select the rigidity accordingly.
147 : class MagneticLens
148 : {
149 :
150 : void updateRigidityBounds(double rigidityMin, double rigidityMax);
151 :
152 : /// Loads part of a lens (one matrix) from file to use it in given rigidity range.
153 : void loadLensPart(const string &filename, double rigidityMin,
154 : double rigidityMax);
155 :
156 : // Stores the individual lenses
157 : std::vector<LensPart*> _lensParts;
158 : Pixelization* _pixelization;
159 : // Checks Matrix, raises Errors if not ok - also generate
160 : // _pixelization if called first time
161 : void _checkMatrix(const ModelMatrixType &M);
162 : // minimum / maximum rigidity that is covered by the lens [Joule]
163 : double _minimumRigidity;
164 : double _maximumRigidity;
165 : static bool _randomSeeded;
166 : double _norm;
167 :
168 : public:
169 : /// Default constructor
170 0 : MagneticLens() :
171 0 : _pixelization(NULL), _minimumRigidity(DBL_MAX), _maximumRigidity(DBL_MIN), _norm(1)
172 : {
173 : }
174 :
175 : /// Constructs lens with predefined healpix order
176 3 : MagneticLens(uint8_t healpixorder) :
177 3 : _pixelization(NULL), _minimumRigidity(DBL_MAX), _maximumRigidity(DBL_MIN)
178 : {
179 3 : _pixelization = new Pixelization(healpixorder);
180 3 : }
181 :
182 : /// Construct lens and load lens from file
183 0 : MagneticLens(const string &filename) :
184 0 : _pixelization(NULL), _minimumRigidity(DBL_MAX), _maximumRigidity(DBL_MIN)
185 : {
186 0 : loadLens(filename);
187 0 : }
188 :
189 : /// Returns the pixelization used
190 : const Pixelization& getPixelization() const
191 : {
192 0 : return (*_pixelization);
193 : }
194 :
195 : /// Default destructor
196 3 : ~MagneticLens()
197 : {
198 3 : if (_pixelization)
199 3 : delete _pixelization;
200 3 : for (std::vector<LensPart*>::iterator iter = _lensParts.begin();
201 6 : iter != _lensParts.end(); iter++)
202 : {
203 3 : delete (*iter);
204 : }
205 : _lensParts.clear();
206 3 : }
207 :
208 : /// Try to transform the comsic ray to a new direction.
209 : /// Returns false and does not change phi and theta if the cosmic ray is
210 : /// lost due to conservation of cosmic ray flux.
211 : /// Rigidity is given in Joule, phi and theta in rad
212 : bool transformCosmicRay(double rigidity, double& phi, double& theta);
213 :
214 : /// Tries transform a cosmic ray with momentum vector p
215 : bool transformCosmicRay(double rigidity, Vector3d &p);
216 :
217 : /// transforms the model array assuming that model points to an array of the
218 : /// correct size. Rigidity is given in Joule
219 : void transformModelVector(double* model, double rigidity) const;
220 :
221 : /// Loads M as part of a lens and use it in given rigidity range with
222 : /// rigidities given in Joule
223 : void setLensPart(const ModelMatrixType &M, double rigidityMin, double rigidityMax);
224 :
225 : /// Loads a lens from a given file, containing lines like
226 : /// lensefile.MLDAT rigidityMin rigidityMax
227 : /// rigidities are given in logarithmic units [log10(E / eV)]
228 : void loadLens(const string &filename);
229 :
230 : /// Normalizes the lens parts to the maximum of sums of columns of
231 : /// every lenspart. By doing this, the lens won't distort the spectrum
232 : void normalizeLens();
233 :
234 : /// Normalizes the lens parts individually. Normalized this way, the
235 : /// lens generally distorts the spectrum of the sources, but deflects
236 : /// the UHECR more efficiently.
237 : void normalizeLensparts();
238 :
239 : /// Checks if rigidity [Joule] is covered by lens
240 : bool rigidityCovered(double rigidity) const;
241 :
242 : /// Normalizes all matrix columns - the lens will then create fake
243 : /// anisotropies, but won't drop particles
244 : void normalizeMatrixColumns();
245 :
246 : /// Returns minimum rigidity covered by lens, in eV
247 : double getMinimumRigidity() const
248 : {
249 0 : return _minimumRigidity / eV;
250 : }
251 : /// Returns maximum rigidity covered by lens, in eV
252 : double getMaximumRigidity() const
253 : {
254 0 : return _maximumRigidity / eV;
255 : }
256 :
257 : // returns the norm used for the lenses
258 : double getNorm()
259 : {
260 0 : return _norm;
261 : }
262 :
263 : /// Returns iterator to the lens part with rigidity Joule
264 : LensPart* getLensPart(double rigidity) const;
265 :
266 : /// Returns all lens parts
267 : const std::vector<LensPart*>& getLensParts() const
268 : {
269 0 : return _lensParts;
270 : }
271 : };
272 :
273 : /** @}*/
274 : } // namespace
275 :
276 : #endif // MAGNETICLENS_HH
|
