Line data Source code
1 : #include "crpropa/InteractionRates.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 :
5 : #include "kiss/path.h"
6 :
7 : #include <vector>
8 : #include <string>
9 : #include <unordered_map>
10 : #include <fstream>
11 :
12 : namespace crpropa {
13 :
14 80 : InteractionRatesHomogeneous::InteractionRatesHomogeneous(std::string rateFile, std::string cumulativeRateFile) {
15 80 : this->ratesName = "interactionRatesHomogeneous";
16 80 : this->isPositionDependent = false;
17 :
18 80 : if (rateFile != "")
19 160 : initRate(rateFile);
20 80 : if (cumulativeRateFile != "")
21 130 : initCumulativeRate(cumulativeRateFile);
22 80 : }
23 :
24 844 : double InteractionRatesHomogeneous::getProcessRate(const double E, const Vector3d &position) const {
25 844 : return interpolate(E, this->tabEnergy, this->tabRate);
26 : }
27 :
28 500 : void InteractionRatesHomogeneous::loadPerformInteractionTabs(const Vector3d &position, std::vector<double> &tabE, std::vector<double> &tabs, std::vector<std::vector<double>> &tabCDF) const {
29 500 : tabE = this->getTabulatedE();
30 500 : tabs = this->getTabulateds();
31 500 : tabCDF = this->getTabulatedCDF();
32 500 : }
33 :
34 80 : void InteractionRatesHomogeneous::setTabulatedEnergy(std::vector<double>& tabEnergy) {
35 80 : this->tabEnergy = tabEnergy;
36 80 : }
37 :
38 80 : void InteractionRatesHomogeneous::setTabulatedRate(std::vector<double>& tabRate) {
39 80 : this->tabRate = tabRate;
40 80 : }
41 :
42 65 : void InteractionRatesHomogeneous::setTabulatedE(std::vector<double>& tabE) {
43 65 : this->tabE = tabE;
44 65 : }
45 :
46 65 : void InteractionRatesHomogeneous::setTabulateds(std::vector<double>& tabs) {
47 65 : this->tabs = tabs;
48 65 : }
49 :
50 65 : void InteractionRatesHomogeneous::setTabulatedCDF(std::vector<std::vector<double>>& tabCDF) {
51 65 : this->tabCDF = tabCDF;
52 65 : }
53 :
54 80 : void InteractionRatesHomogeneous::initRate(std::string filename){
55 80 : if(is_directory(filename))
56 0 : throw std::runtime_error("InteractionRatesHomogeneous: The given filename " + filename + " is a directory!\
57 0 : If you wanted to use position dependent photon fields instead, set the correct InteractionRates class first!");
58 :
59 80 : std::ifstream infile(filename.c_str());
60 80 : if (!infile.good())
61 0 : throw std::runtime_error("InteractionRatesHomogeneous: could not open file " + filename);
62 :
63 : std::vector<double> tabEnergy;
64 : std::vector<double> tabRate;
65 :
66 19035 : while (infile.good()) {
67 :
68 18955 : if (infile.peek() != '#') {
69 : double a, b;
70 : infile >> a >> b;
71 18635 : if (infile) {
72 18555 : tabEnergy.push_back(pow(10, a) * eV);
73 18555 : tabRate.push_back(b / Mpc);
74 : }
75 : }
76 18955 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
77 :
78 : }
79 :
80 80 : infile.close();
81 :
82 80 : setTabulatedEnergy(tabEnergy);
83 80 : setTabulatedRate(tabRate);
84 80 : }
85 :
86 65 : void InteractionRatesHomogeneous::initCumulativeRate(std::string filename){
87 65 : if(is_directory(filename))
88 0 : throw std::runtime_error("InteractionRatesHomogeneous: The given filename " + filename + " is a directory!\
89 0 : If you wanted to use position dependent photon fields instead, set the correct InteractionRates class first!");
90 :
91 65 : std::ifstream infile(filename.c_str());
92 :
93 65 : if (!infile.good())
94 0 : throw std::runtime_error("InteractionRatesHomogeneous: could not open file " + filename);
95 :
96 : std::vector<double> tabE;
97 : std::vector<double> tabs;
98 : std::vector<std::vector<double>> tabCDF;
99 :
100 : // skip header
101 325 : while (infile.peek() == '#')
102 260 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
103 :
104 : // read s values in first line
105 : double a;
106 : infile >> a; // skip first value
107 8253 : while (infile.good() and (infile.peek() != '\n')) {
108 : infile >> a;
109 8188 : tabs.push_back(pow(10, a) * eV * eV);
110 : }
111 :
112 : // read all following lines: E, cdf values
113 15404 : while (infile.good()) {
114 : infile >> a;
115 15404 : if (!infile)
116 : break; // end of file
117 15339 : tabE.push_back(pow(10, a) * eV);
118 : std::vector<double> cdf;
119 1999132 : for (int i = 0; i < tabs.size(); i++) {
120 : infile >> a;
121 1983793 : cdf.push_back(a / Mpc);
122 : }
123 15339 : tabCDF.push_back(cdf);
124 15339 : }
125 65 : infile.close();
126 :
127 65 : setTabulatedE(tabE);
128 65 : setTabulateds(tabs);
129 65 : setTabulatedCDF(tabCDF);
130 65 : }
131 :
132 :
133 0 : InteractionRatesPositionDependent::InteractionRatesPositionDependent(
134 0 : std::string rateFilePath, std::string cumulativeRateFilePath, ref_ptr<Surface> surface) {
135 :
136 0 : this->ratesName = "interactionRatesPositionDependent";
137 0 : this->isPositionDependent = true;
138 0 : this->surface = surface;
139 :
140 0 : if (rateFilePath != "")
141 0 : initRate(rateFilePath);
142 0 : if (cumulativeRateFilePath != "")
143 0 : initCumulativeRate(cumulativeRateFilePath);
144 0 : }
145 :
146 0 : int InteractionRatesPositionDependent::findClosestGridPoint(const Vector3d &position) const {
147 0 : if (!tree) {
148 0 : throw std::runtime_error("KD-Tree not initialized!");
149 : }
150 :
151 : unsigned int closestIndex;
152 : double closestDistSquared;
153 0 : double queryPoint[3] = { position.x, position.y, position.z };
154 :
155 : this->tree->knnSearch(queryPoint, 1, &closestIndex, &closestDistSquared);
156 0 : return this->cloud.ids[closestIndex];
157 :
158 : }
159 :
160 0 : std::vector<double> InteractionRatesPositionDependent::getClosestRate(const Vector3d &position) const {
161 0 : int iMin = findClosestGridPoint(position);
162 0 : return tabRate[iMin];
163 : }
164 :
165 0 : std::vector<double> InteractionRatesPositionDependent::getClosests(const Vector3d &position) const {
166 0 : int iMin = findClosestGridPoint(position);
167 0 : return tabs[iMin];
168 : }
169 :
170 0 : std::vector<std::vector<double>> InteractionRatesPositionDependent::getClosestCDF(const Vector3d &position) const {
171 0 : int iMin = findClosestGridPoint(position);
172 0 : return tabCDF[iMin];
173 : }
174 :
175 0 : double InteractionRatesPositionDependent::getProcessRate(const double E, const Vector3d &position) const {
176 0 : std::vector<double> tabRate = this->getClosestRate(position);
177 : // compute the interaction rate for the given candidate energy, E
178 0 : return interpolate(E, this->tabEnergy, tabRate);
179 0 : }
180 :
181 0 : void InteractionRatesPositionDependent::loadPerformInteractionTabs(const Vector3d &position, std::vector<double> &tabE, std::vector<double> &tabs, std::vector<std::vector<double>> &tabCDF) const {
182 :
183 : std::vector<double> E = this->getTabulatedE();
184 0 : std::vector<double> s = this->getClosests(position);
185 0 : std::vector<std::vector<double>> CDF = this->getClosestCDF(position);
186 :
187 0 : tabE = E;
188 0 : tabs = s;
189 0 : tabCDF = CDF;
190 0 : }
191 :
192 0 : void InteractionRatesPositionDependent::setTabulatedEnergy(std::vector<double>& tabEnergy) {
193 0 : this->tabEnergy = tabEnergy;
194 0 : }
195 :
196 0 : void InteractionRatesPositionDependent::setTabulatedRate(std::vector<std::vector<double>>& tabRate) {
197 0 : this->tabRate = tabRate;
198 0 : }
199 :
200 0 : void InteractionRatesPositionDependent::setTabulatedE(std::vector<double>& tabE) {
201 0 : this->tabE = tabE;
202 0 : }
203 :
204 0 : void InteractionRatesPositionDependent::setTabulateds(std::vector<std::vector<double>>& tabs) {
205 0 : this->tabs = tabs;
206 0 : }
207 :
208 0 : void InteractionRatesPositionDependent::setTabulatedCDF(std::vector<std::vector<std::vector<double>>>& tabCDF) {
209 0 : this->tabCDF = tabCDF;
210 0 : }
211 :
212 0 : void InteractionRatesPositionDependent::setPhotonDict(std::unordered_map<int, Vector3d>& photonDict) {
213 : this->photonDict = photonDict;
214 :
215 : // delete old clouds
216 : this->cloud.points.clear();
217 : this->cloud.ids.clear();
218 :
219 0 : for (const auto& el : this->photonDict) {
220 0 : this->cloud.ids.push_back(el.first);
221 0 : this->cloud.points.push_back(el.second);
222 : }
223 :
224 : // delete old tree
225 0 : if (this->tree) {
226 0 : delete this->tree;
227 : }
228 :
229 : int maxLeafTree = 20;
230 : int nThreads = 4;
231 : nanoflann::KDTreeSingleIndexAdaptorFlags flag;
232 :
233 0 : this->tree = new KDTree(3, cloud, nanoflann::KDTreeSingleIndexAdaptorParams(maxLeafTree, flag, nThreads));
234 0 : this->tree->buildIndex();
235 0 : }
236 :
237 0 : void InteractionRatesPositionDependent::setSurface(ref_ptr<Surface> surface) {
238 0 : this->surface = surface;
239 0 : }
240 :
241 0 : void InteractionRatesPositionDependent::initRate(std::string filepath){
242 0 : if(!is_directory(filepath))
243 0 : throw std::runtime_error("InteractionRatesPositionDependent: The given path " + filepath + " is not a directory or does not exist!\
244 0 : If you wanted to use homogeneous photon fields instead, set the correct InteractionRates class first!");
245 :
246 : std::vector<std::vector<double>> tabRate;
247 :
248 : std::unordered_map<int, Vector3d> photonDict;
249 0 : int iFile = 0;
250 :
251 : std::vector<std::string> dirs;
252 0 : if(!list_directory(filepath, dirs))
253 0 : throw std::runtime_error("Could not find any files in " + filepath + "!\n");
254 :
255 0 : for (auto const& filename : dirs) {
256 :
257 : // the input filename here should be a string
258 : //check if it is correct, i.e. a proper filename string
259 0 : std::ifstream infile(concat_path(filepath, filename).c_str());
260 0 : if (!infile.good())
261 0 : throw std::runtime_error("InteractionRatesPositionDependent: could not open file " + concat_path(filepath, filename));
262 :
263 : std::vector<double> vecEnergy;
264 : std::vector<double> vecRate;
265 :
266 : double x, y, z;
267 : std::string str;
268 0 : std::stringstream ss;
269 :
270 0 : std::string filename_split = splitFilename(filename);
271 : ss << filename_split;
272 :
273 : int iLine = 0;
274 :
275 0 : std::locale::global(std::locale("C"));
276 :
277 0 : while (getline(ss, str, '_')) {
278 0 : if (iLine == 3) {
279 0 : x = std::stod(str) * kpc;
280 : }
281 0 : if (iLine == 4) {
282 0 : y = std::stod(str) * kpc;
283 : }
284 0 : if (iLine == 5) {
285 0 : z = std::stod(str) * kpc;
286 : }
287 0 : iLine = iLine + 1;
288 : }
289 :
290 : Vector3d vPos(x, y, z);
291 :
292 : // continue when not "inside" surface
293 0 : if (getSurface() && getSurface()->distance(vPos)>=0)
294 : continue;
295 :
296 : photonDict[iFile] = vPos;
297 :
298 0 : while (infile.good()) {
299 0 : if (infile.peek() != '#') {
300 : double a, b;
301 : infile >> a >> b;
302 0 : if (infile) {
303 0 : if (iFile == 0) {
304 0 : vecEnergy.push_back(pow(10, a) * eV);
305 0 : setTabulatedEnergy(vecEnergy);
306 : }
307 0 : vecRate.push_back(b / Mpc);
308 : }
309 : }
310 0 : infile.ignore(std::numeric_limits< std::streamsize >::max(), '\n');
311 : }
312 :
313 0 : tabRate.push_back(vecRate);
314 :
315 0 : iFile = iFile + 1;
316 0 : infile.close();
317 0 : }
318 :
319 0 : if (tabRate.empty())
320 0 : throw std::runtime_error("Rate's table empty! Check if the surface is properly set.");
321 :
322 0 : setTabulatedRate(tabRate);
323 0 : setPhotonDict(photonDict);
324 0 : }
325 :
326 0 : void InteractionRatesPositionDependent::initCumulativeRate(std::string filepath){
327 0 : if(!is_directory(filepath))
328 0 : throw std::runtime_error("InteractionRatesPositionDependent: The given path " + filepath + " is not a directory or does not exist!\
329 0 : If you wanted to use homogeneous photon fields instead, set the correct InteractionRates class first!");
330 :
331 : std::vector<std::vector<double>> tabs;
332 : std::vector<std::vector<std::vector<double>>> tabCDF;
333 :
334 : int iFile = 0;
335 : std::vector<std::string> dirs;
336 0 : if(!list_directory(filepath, dirs))
337 0 : throw std::runtime_error("Could not find any files in " + filepath + "!\n");
338 :
339 0 : for (auto const& filename : dirs) {
340 :
341 : std::vector<double> vecE;
342 : std::vector<double> vecs;
343 : std::vector<std::vector<double>> vecCDF;
344 :
345 0 : std::ifstream infile(concat_path(filepath, filename).c_str());
346 :
347 0 : if (!infile.good())
348 0 : throw std::runtime_error("InteractionRatesPositionDependent: could not open file " + concat_path(filepath, filename));
349 :
350 : double x, y, z;
351 : std::string str;
352 0 : std::stringstream ss;
353 :
354 0 : std::string filename_split = splitFilename(filename);
355 : ss << filename_split;
356 :
357 : int iLine = 0;
358 :
359 0 : std::locale::global(std::locale("C"));
360 :
361 0 : while (getline(ss, str, '_')) {
362 0 : if (iLine == 3) {
363 0 : x = -std::stod(str) * kpc;
364 : }
365 0 : if (iLine == 4) {
366 0 : y = std::stod(str) * kpc;
367 : }
368 0 : if (iLine == 5) {
369 0 : z = std::stod(str) * kpc;
370 : }
371 0 : iLine = iLine + 1;
372 : }
373 :
374 : Vector3d vPos(x, y, z);
375 :
376 : // continue when not "inside" surface
377 0 : if (getSurface() && getSurface()->distance(vPos)>=0)
378 : continue;
379 :
380 : // skip header
381 0 : while (infile.peek() == '#')
382 0 : infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
383 :
384 : // read s values in first line
385 : double a;
386 : infile >> a; // skip first value
387 0 : while (infile.good() and (infile.peek() != '\n')) {
388 : infile >> a;
389 0 : vecs.push_back(pow(10, a) * eV * eV);
390 : }
391 :
392 : // read all following lines: E, cdf values
393 0 : while (infile.good()) {
394 : infile >> a;
395 0 : if (!infile)
396 : break; // end of file
397 0 : if (iFile == 0) {
398 0 : vecE.push_back(pow(10, a) * eV);
399 0 : setTabulatedE(vecE);
400 : }
401 : std::vector<double> cdf;
402 0 : for (int i = 0; i < tabs.size(); i++) {
403 : infile >> a;
404 0 : cdf.push_back(a / Mpc);
405 : }
406 0 : vecCDF.push_back(cdf);
407 0 : }
408 :
409 0 : iFile = iFile + 1;
410 :
411 0 : tabs.push_back(vecs);
412 0 : tabCDF.push_back(vecCDF);
413 0 : infile.close();
414 0 : }
415 :
416 0 : setTabulateds(tabs);
417 0 : setTabulatedCDF(tabCDF);
418 0 : }
419 :
420 : } //namespace crpropa
|
