Line data Source code
1 : #ifndef CRPROPA_INTERACTIONRATES_H
2 : #define CRPROPA_INTERACTIONRATES_H
3 :
4 : #include "crpropa/Common.h"
5 : #include "crpropa/Referenced.h"
6 : #include "crpropa/Vector3.h"
7 : #include "crpropa/Geometry.h"
8 :
9 : #include "nanoflann/nanoflann.hpp"
10 :
11 : #include <vector>
12 : #include <string>
13 : #include <unordered_map>
14 :
15 : namespace crpropa {
16 :
17 0 : struct PointCloud {
18 : std::vector<Vector3d> points;
19 : std::vector<int> ids;
20 :
21 : inline size_t kdtree_get_point_count() const { return points.size(); }
22 :
23 : inline double kdtree_get_pt(const size_t idx, const size_t dim) const {
24 0 : if (dim == 0) return points[idx].x;
25 0 : if (dim == 1) return points[idx].y;
26 0 : return points[idx].z;
27 : }
28 :
29 : // optional bounding-box computation (required by nanoflann)
30 : template <class BBOX>
31 : bool kdtree_get_bbox(BBOX& /*bb*/) const {
32 : return false; // no bounding box optimization
33 : }
34 :
35 : };
36 :
37 : using KDTree = nanoflann::KDTreeSingleIndexAdaptor<
38 : nanoflann::L2_Simple_Adaptor<double, PointCloud>,
39 : PointCloud,
40 : 3
41 : >;
42 :
43 : /**
44 : * \addtogroup InteractionRates
45 : * @{
46 : */
47 :
48 : /**
49 : @class Interaction Rates
50 : @brief Abstract base class for photon fields interaction rates.
51 : */
52 80 : class InteractionRates: public Referenced {
53 : public:
54 : virtual double getProcessRate(const double E, const Vector3d &position) const = 0;
55 : virtual void loadPerformInteractionTabs(
56 : const Vector3d &position,
57 : std::vector<double> &tabE,
58 : std::vector<double> &tabs,
59 : std::vector<std::vector<double>> &tabCDF
60 : ) const = 0;
61 :
62 : std::string getRatesName() const {
63 0 : return this->ratesName;
64 : }
65 :
66 : bool hasPositionDependence() const {
67 0 : return this->isPositionDependent;
68 : }
69 :
70 : void setRatesName(std::string ratesName) {
71 0 : this->ratesName = ratesName;
72 0 : }
73 :
74 : virtual void initRate(std::string path) = 0;
75 : virtual void initCumulativeRate(std::string path) = 0;
76 :
77 : protected:
78 :
79 : std::string ratesName = "AbstractInteractionRates";
80 : bool isPositionDependent = false;
81 :
82 : };
83 :
84 : /**
85 : @class InteractionRateHomogeneous
86 : @brief Interaction rates decorator for tabulated homogeneous photon fields.
87 :
88 : This class handles the interaction rates for homogeneous photon fields.
89 : It is mainly used in the EM* interaction classes.
90 : For position dependend photon fields use InteractionRatesPositionDependent.
91 : */
92 : class InteractionRatesHomogeneous: public InteractionRates {
93 : public:
94 : /** Constructor of InteractionRatesHomogeneous
95 : * @param RateFile Path to the file containing the interaction rate data
96 : * @param CumulativeRateFile Path to the file containing the cumulative interaction rate data
97 : */
98 : InteractionRatesHomogeneous(std::string RateFile = "", std::string CumulativeRateFile = "");
99 :
100 0 : inline std::vector<double> getTabulatedEnergy() const {return tabEnergy;}
101 0 : inline std::vector<double> getTabulatedRate() const {return tabRate;}
102 500 : inline std::vector<double> getTabulatedE() const {return tabE;}
103 500 : inline std::vector<double> getTabulateds() const {return tabs;}
104 500 : inline std::vector<std::vector<double>> getTabulatedCDF() const {return tabCDF;}
105 :
106 : /** Returns the interaction rate for the given candidate energy E
107 : * this is just a wrapper for interpolate(E, tabEnergy, tabRate) and is not positon dependend.
108 : * @param E Candidate Energy
109 : * @param position Position is never used here
110 : */
111 : double getProcessRate(const double E, const Vector3d &position = Vector3d(0,0,0)) const;
112 : /** Loads tabE, tabs and tabCDF
113 : * Alternatively you could use getTabulated*() functions directly here since nothing is position dependend here
114 : * @param position Position is not used here
115 : * @param tabE Vector to save tabulated energies to
116 : * @param tabs Vector to save tabulated rates to
117 : * @param tabCDF Vector to save tabulated CDF to
118 : */
119 : void loadPerformInteractionTabs(
120 : const Vector3d &position,
121 : std::vector<double> &tabE,
122 : std::vector<double> &tabs,
123 : std::vector<std::vector<double>> &tabCDF
124 : ) const;
125 :
126 : void setTabulatedEnergy (std::vector<double>& tabEnergy);
127 : void setTabulatedRate (std::vector<double>& tabRate);
128 : void setTabulatedE (std::vector<double>& tabE);
129 : void setTabulateds (std::vector<double>& tabs);
130 : void setTabulatedCDF (std::vector<std::vector<double>>& tabCDF);
131 :
132 : /** Loads the interaction rate
133 : * This function loads the interaction rate from a given filename
134 : * @param filename The name of the file containing the interaction rates
135 : */
136 : void initRate(std::string filename);
137 : /** Loads the cumulative interaction rate
138 : * This function loads the interaction rate from a given filename
139 : * @param filename The name of the file containing the interaction rates
140 : */
141 : void initCumulativeRate(std::string filename);
142 :
143 : protected:
144 :
145 : // tabulated interaction rates 1/lambda(E)
146 : std::vector<double> tabEnergy; //!< electron energy in [J]
147 : std::vector<double> tabRate; //!< interaction rate in [1/m]
148 :
149 : // tabulated CDF(s_kin, E) = cumulative differential interaction rate
150 : std::vector<double> tabE; //!< electron energy in [J]
151 : std::vector<double> tabs; //!< s_kin = s - m^2 in [J**2]
152 : std::vector<std::vector<double>> tabCDF; //!< cumulative interaction rate
153 :
154 : };
155 :
156 : /**
157 : @class InteractionRatePositionDependent
158 : @brief Interaction rates decorator for tabulated position dependent photon fields.
159 : */
160 : class InteractionRatesPositionDependent: public InteractionRates {
161 : private:
162 : /** Finds closest point cloud ID
163 : * This function takes a position vector and searches for the closes Grid point,
164 : * it returns the point cloud ID to be used
165 : * @param position The position to search for a grid point
166 : */
167 : int findClosestGridPoint(const Vector3d &position) const;
168 :
169 : public:
170 : /** Constructor of InteractionRatesPositionDependent
171 : * @param RateFilePath Path containing the interaction rates files (* /Rate)
172 : * @param CumulativeRateFilePath Path containing the cumulative interaction rates files (* /CumulativeRate)
173 : * @param surface Closed surface to confine the grid nodes to be uploaded (optional)
174 : */
175 : InteractionRatesPositionDependent(
176 : std::string RateFilePath = "",
177 : std::string CumulativeRateFilePath = "",
178 : ref_ptr<Surface> surface = NULL
179 : );
180 :
181 0 : inline std::vector<double> getTabulatedEnergy() const {return tabEnergy;}
182 0 : inline std::vector<std::vector<double>> getTabulatedRate() const {return tabRate;}
183 0 : inline std::vector<double> getTabulatedE() const {return tabE;}
184 0 : inline std::vector<std::vector<double>> getTabulateds() const {return tabs;}
185 0 : inline std::vector<std::vector<std::vector<double>>> getTabulatedCDF() const {return tabCDF;}
186 : inline std::unordered_map<int, Vector3d> getPhotonDict() const {return photonDict;}
187 : std::vector<double> getClosestRate(const Vector3d &position) const;
188 : std::vector<double> getClosests(const Vector3d &position) const;
189 : std::vector<std::vector<double>> getClosestCDF(const Vector3d &position) const;
190 :
191 : /** Returns the interaction rate for the given candidate energy E and position
192 : * This function calculated the interaction rate dependend on the given postion and particle energy E
193 : * @param E Candidate Energy
194 : * @param position Current position of particle
195 : */
196 : double getProcessRate(const double E, const Vector3d &position) const;
197 : /** Loads tabE, tabs and tabCDF dependend on the position
198 : * Loads the tabulated Energies, rates and CDF dependend on the given position
199 : * @param position Position to load tab* for
200 : * @param tabE Vector to save tabulated energies to
201 : * @param tabs Vector to save tabulated rates to
202 : * @param tabCDF Vector to save tabulated CDF to
203 : */
204 : void loadPerformInteractionTabs(
205 : const Vector3d &position,
206 : std::vector<double> &tabE,
207 : std::vector<double> &tabs,
208 : std::vector<std::vector<double>> &tabCDF
209 : ) const;
210 :
211 : void setTabulatedEnergy (std::vector<double>& tabEnergy);
212 : void setTabulatedRate (std::vector<std::vector<double>>& tabRate);
213 : void setTabulatedE (std::vector<double>& tabE);
214 : void setTabulateds (std::vector<std::vector<double>>& tabs);
215 : void setTabulatedCDF (std::vector<std::vector<std::vector<double>>>& tabCDF);
216 : void setPhotonDict (std::unordered_map<int, Vector3d>& photonDict);
217 :
218 : /** Apply a surface that confine the position dependent photon field
219 : * @param surface Closed surface to confine the grid nodes to be uploaded
220 : */
221 : void setSurface(ref_ptr<Surface> surface);
222 : inline ref_ptr<Surface> getSurface() const {return surface;}
223 :
224 : /** Loads the interaction rate
225 : * This function loads the position dependent interaction rate
226 : * @param filepath The name of the folder containing the interaction rates (* /Rate)
227 : */
228 : void initRate(std::string filepath);
229 : /** Loads the interaction rate
230 : * This function loads the cumulative position dependent interaction rate
231 : * @param filepath The name of the folder containing the cumulative interaction rates (* /CumulativeRate)
232 : */
233 : void initCumulativeRate(std::string filepath);
234 :
235 : protected:
236 :
237 : // tabulated interaction rates 1/lambda(E)
238 : std::vector<double> tabEnergy; //!< electron energy in [J], assuming the same energy binning in each node
239 : std::vector<std::vector<double>> tabRate; //!< interaction rate in [1/m]
240 :
241 : // tabulated CDF(s_kin, E) = cumulative differential interaction rate
242 : std::vector<double> tabE; //!< electron energy in [J], assuming the same energy binning in each node
243 : std::vector<std::vector<double>> tabs; //!< s_kin = s - m^2 in [J**2]
244 : std::vector<std::vector<std::vector<double>>> tabCDF; //!< cumulative interaction rate
245 : std::unordered_map<int, Vector3d> photonDict; //!< dictionary to link tables to spatial coordinates
246 :
247 : PointCloud cloud; //!< point cloud for nanoflann KD-tree
248 : KDTree* tree = nullptr; //!< pointer to the KD Tree
249 : ref_ptr<Surface> surface;
250 :
251 : };
252 :
253 : } // namespace crpropa
254 :
255 : #endif // CRPROPA_INTERACTIONRATES_H
|
