Line data Source code
1 : #include "crpropa/EmissionMap.h"
2 : #include "crpropa/Random.h"
3 : #include "crpropa/Units.h"
4 :
5 : #include "kiss/logger.h"
6 :
7 : #include <fstream>
8 :
9 : namespace crpropa {
10 :
11 0 : CylindricalProjectionMap::CylindricalProjectionMap() : nPhi(360), nTheta(180), dirty(false), pdf(nPhi* nTheta, 0), cdf(nPhi* nTheta, 0) {
12 0 : sPhi = 2. * M_PI / nPhi;
13 0 : sTheta = 2. / nTheta;
14 0 : }
15 :
16 6 : CylindricalProjectionMap::CylindricalProjectionMap(size_t nPhi, size_t nTheta) : nPhi(nPhi), nTheta(nTheta), dirty(false), pdf(nPhi* nTheta, 0), cdf(nPhi* nTheta, 0) {
17 6 : sPhi = 2 * M_PI / nPhi;
18 6 : sTheta = 2. / nTheta;
19 6 : }
20 :
21 4 : void CylindricalProjectionMap::fillBin(const Vector3d& direction, double weight) {
22 4 : size_t bin = binFromDirection(direction);
23 4 : fillBin(bin, weight);
24 4 : }
25 :
26 129604 : void CylindricalProjectionMap::fillBin(size_t bin, double weight) {
27 129604 : pdf[bin] += weight;
28 129604 : dirty = true;
29 129604 : }
30 :
31 1 : Vector3d CylindricalProjectionMap::drawDirection() const {
32 1 : if (dirty)
33 1 : updateCdf();
34 :
35 1 : size_t bin = Random::instance().randBin(cdf);
36 :
37 1 : return directionFromBin(bin);
38 : }
39 :
40 0 : bool CylindricalProjectionMap::checkDirection(const Vector3d &direction) const {
41 0 : size_t bin = binFromDirection(direction);
42 0 : return pdf[bin];
43 : }
44 :
45 :
46 0 : const std::vector<double>& CylindricalProjectionMap::getPdf() const {
47 0 : return pdf;
48 : }
49 :
50 5 : std::vector<double>& CylindricalProjectionMap::getPdf() {
51 5 : return pdf;
52 : }
53 :
54 0 : const std::vector<double>& CylindricalProjectionMap::getCdf() const {
55 0 : return cdf;
56 : }
57 :
58 0 : size_t CylindricalProjectionMap::getNPhi() {
59 0 : return nPhi;
60 : }
61 :
62 0 : size_t CylindricalProjectionMap::getNTheta() {
63 0 : return nTheta;
64 : }
65 :
66 : /*
67 : * Cylindrical Coordinates
68 : * iPhi -> [0, 2*pi]
69 : * iTheta -> [0, 2]
70 : *
71 : * Spherical Coordinates
72 : * phi -> [-pi, pi]
73 : * theta -> [0, pi]
74 : */
75 6 : size_t CylindricalProjectionMap::binFromDirection(const Vector3d& direction) const {
76 : // convert to cylindrical
77 6 : double phi = direction.getPhi() + M_PI;
78 6 : double theta = sin(M_PI_2 - direction.getTheta()) + 1;
79 :
80 : // to indices
81 6 : size_t iPhi = phi / sPhi;
82 6 : size_t iTheta = theta / sTheta;
83 :
84 : // interleave
85 6 : size_t bin = iTheta * nPhi + iPhi;
86 6 : return bin;
87 : }
88 :
89 2 : Vector3d CylindricalProjectionMap::directionFromBin(size_t bin) const {
90 : // deinterleave
91 2 : double iPhi = bin % nPhi;
92 2 : double iTheta = bin / nPhi;
93 :
94 : // any where in the bin
95 2 : iPhi += Random::instance().rand();
96 2 : iTheta += Random::instance().rand();
97 :
98 : // cylindrical Coordinates
99 2 : double phi = iPhi * sPhi;
100 2 : double theta = iTheta * sTheta;
101 :
102 : // sphericala Coordinates
103 2 : phi = phi - M_PI;
104 2 : theta = M_PI_2 - asin(theta - 1);
105 :
106 : Vector3d v;
107 2 : v.setRThetaPhi(1.0, theta, phi);
108 2 : return v;
109 : }
110 :
111 1 : void CylindricalProjectionMap::updateCdf() const {
112 1 : if (dirty) {
113 1 : cdf[0] = pdf[0];
114 64800 : for (size_t i = 1; i < pdf.size(); i++) {
115 64799 : cdf[i] = cdf[i-1] + pdf[i];
116 : }
117 1 : dirty = false;
118 : }
119 1 : }
120 :
121 2 : EmissionMap::EmissionMap() : minEnergy(0.0001 * EeV), maxEnergy(10000 * EeV),
122 2 : nEnergy(8*2), nPhi(360), nTheta(180) {
123 2 : logStep = log10(maxEnergy / minEnergy) / nEnergy;
124 2 : }
125 :
126 0 : EmissionMap::EmissionMap(size_t nPhi, size_t nTheta, size_t nEnergy) : minEnergy(0.0001 * EeV), maxEnergy(10000 * EeV),
127 0 : nEnergy(nEnergy), nPhi(nPhi), nTheta(nTheta) {
128 0 : logStep = log10(maxEnergy / minEnergy) / nEnergy;
129 0 : }
130 :
131 1 : EmissionMap::EmissionMap(size_t nPhi, size_t nTheta, size_t nEnergy, double minEnergy, double maxEnergy) : minEnergy(minEnergy), maxEnergy(maxEnergy), nEnergy(nEnergy), nPhi(nPhi), nTheta(nTheta) {
132 1 : logStep = log10(maxEnergy / minEnergy) / nEnergy;
133 1 : }
134 :
135 1 : double EmissionMap::energyFromBin(size_t bin) const {
136 1 : return pow(10, log10(minEnergy) + logStep * bin);
137 : }
138 :
139 11 : size_t EmissionMap::binFromEnergy(double energy) const {
140 11 : return log10(energy / minEnergy) / logStep;
141 : }
142 :
143 4 : void EmissionMap::fillMap(int pid, double energy, const Vector3d& direction, double weight) {
144 4 : getMap(pid, energy)->fillBin(direction, weight);
145 4 : }
146 :
147 0 : void EmissionMap::fillMap(const ParticleState& state, double weight) {
148 0 : fillMap(state.getId(), state.getEnergy(), state.getDirection(), weight);
149 0 : }
150 :
151 1 : EmissionMap::map_t &EmissionMap::getMaps() {
152 1 : return maps;
153 : }
154 :
155 2 : const EmissionMap::map_t &EmissionMap::getMaps() const {
156 2 : return maps;
157 : }
158 :
159 3 : bool EmissionMap::drawDirection(int pid, double energy, Vector3d& direction) const {
160 3 : key_t key(pid, binFromEnergy(energy));
161 : map_t::const_iterator i = maps.find(key);
162 :
163 3 : if (i == maps.end() || !i->second.valid()) {
164 : return false;
165 : } else {
166 1 : direction = i->second->drawDirection();
167 1 : return true;
168 : }
169 : }
170 :
171 0 : bool EmissionMap::drawDirection(const ParticleState& state, Vector3d& direction) const {
172 0 : return drawDirection(state.getId(), state.getEnergy(), direction);
173 : }
174 :
175 0 : bool EmissionMap::checkDirection(int pid, double energy, const Vector3d& direction) const {
176 0 : key_t key(pid, binFromEnergy(energy));
177 : map_t::const_iterator i = maps.find(key);
178 :
179 0 : if (i == maps.end() || !i->second.valid()) {
180 : return false;
181 : } else {
182 0 : return i->second->checkDirection(direction);
183 : }
184 : }
185 :
186 0 : bool EmissionMap::checkDirection(const ParticleState& state) const {
187 0 : return checkDirection(state.getId(), state.getEnergy(), state.getDirection());
188 : }
189 :
190 0 : bool EmissionMap::hasMap(int pid, double energy) {
191 0 : key_t key(pid, binFromEnergy(energy));
192 : map_t::iterator i = maps.find(key);
193 0 : if (i == maps.end() || !i->second.valid())
194 : return false;
195 : else
196 0 : return true;
197 : }
198 :
199 7 : ref_ptr<CylindricalProjectionMap> EmissionMap::getMap(int pid, double energy) {
200 7 : key_t key(pid, binFromEnergy(energy));
201 : map_t::iterator i = maps.find(key);
202 7 : if (i == maps.end() || !i->second.valid()) {
203 5 : ref_ptr<CylindricalProjectionMap> cpm = new CylindricalProjectionMap(nPhi, nTheta);
204 5 : maps[key] = cpm;
205 : return cpm;
206 : } else {
207 : return i->second;
208 : }
209 : }
210 :
211 0 : void EmissionMap::save(const std::string &filename) {
212 0 : std::ofstream out(filename.c_str());
213 0 : out.imbue(std::locale("C"));
214 :
215 0 : for (map_t::iterator i = maps.begin(); i != maps.end(); i++) {
216 0 : if (!i->second.valid())
217 0 : continue;
218 0 : out << i->first.first << " " << i->first.second << " " << energyFromBin(i->first.second) << " ";
219 0 : out << i->second->getNPhi() << " " << i->second->getNTheta();
220 0 : const std::vector<double> &pdf = i->second->getPdf();
221 0 : for (size_t i = 0; i < pdf.size(); i++)
222 0 : out << " " << pdf[i];
223 : out << std::endl;
224 : }
225 0 : }
226 :
227 1 : void EmissionMap::merge(const EmissionMap *other) {
228 1 : if (other == 0)
229 : return;
230 1 : map_t::const_iterator i = other->getMaps().begin();
231 1 : map_t::const_iterator end = other->getMaps().end();
232 3 : for(;i != end; i++) {
233 2 : if (!i->second.valid())
234 0 : continue;
235 :
236 2 : std::vector<double> &otherpdf = i->second->getPdf();
237 2 : ref_ptr<CylindricalProjectionMap> cpm = getMap(i->first.first, i->first.second);
238 :
239 2 : if (otherpdf.size() != cpm->getPdf().size()) {
240 0 : throw std::runtime_error("PDF size mismatch!");
241 : break;
242 : }
243 :
244 129602 : for (size_t k = 0; k < otherpdf.size(); k++) {
245 129600 : cpm->fillBin(k, otherpdf[k]);
246 : }
247 : }
248 : }
249 :
250 0 : void EmissionMap::merge(const std::string &filename) {
251 0 : EmissionMap em;
252 0 : em.load(filename);
253 0 : merge(&em);
254 0 : }
255 :
256 0 : void EmissionMap::load(const std::string &filename) {
257 0 : std::ifstream in(filename.c_str());
258 0 : in.imbue(std::locale("C"));
259 :
260 0 : while(in.good()) {
261 0 : key_t key;
262 : double tmp;
263 : size_t nPhi_, nTheta_;
264 0 : in >> key.first >> key.second >> tmp;
265 : in >> nPhi_ >> nTheta_;
266 :
267 0 : if (!in.good()) {
268 0 : KISS_LOG_ERROR << "Invalid line: " << key.first << " " << key.second << " " << nPhi_ << " " << nTheta_;
269 0 : break;
270 : }
271 :
272 0 : if (nPhi != nPhi_) {
273 0 : KISS_LOG_WARNING << "nPhi mismatch: " << nPhi << " " << nPhi_;
274 : }
275 0 : if (nTheta != nTheta_) {
276 0 : KISS_LOG_WARNING << "nTheta mismatch: " << nTheta << " " << nTheta_;
277 : }
278 :
279 0 : ref_ptr<CylindricalProjectionMap> cpm = new CylindricalProjectionMap(nPhi_, nTheta_);
280 0 : std::vector<double> &pdf = cpm->getPdf();
281 0 : for (size_t i = 0; i < pdf.size(); i++)
282 : in >> pdf[i];
283 :
284 0 : if (in.good()) {
285 0 : maps[key] = cpm;
286 : // std::cout << "added " << key.first << " " << key.second << std::endl;
287 : } else {
288 0 : KISS_LOG_WARNING << "Invalid data in line: " << key.first << " " << key.second << " " << nPhi_ << " " << nTheta_ << "\n";
289 : }
290 : }
291 :
292 0 : }
293 :
294 : } // namespace crpropa
|
