Line data Source code
1 : #include "crpropa/module/Observer.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/ParticleID.h"
4 : #include "crpropa/Cosmology.h"
5 :
6 : #include "kiss/logger.h"
7 :
8 : #include <iostream>
9 : #include <cmath>
10 :
11 : namespace crpropa {
12 :
13 : // Observer -------------------------------------------------------------------
14 10 : Observer::Observer() :
15 10 : makeInactive(true), clone(false) {
16 10 : }
17 :
18 10 : void Observer::add(ObserverFeature *feature) {
19 10 : features.push_back(feature);
20 10 : }
21 :
22 2 : void Observer::onDetection(Module *action, bool clone_) {
23 2 : detectionAction = action;
24 2 : clone = clone_;
25 2 : }
26 :
27 476 : void Observer::process(Candidate *candidate) const {
28 : // loop over all features and have them check the particle
29 : DetectionState state = NOTHING;
30 952 : for (int i = 0; i < features.size(); i++) {
31 476 : DetectionState s = features[i]->checkDetection(candidate);
32 476 : if (s == VETO)
33 : state = VETO;
34 476 : else if ((s == DETECTED) && (state != VETO))
35 : state = DETECTED;
36 : }
37 :
38 476 : if (state == DETECTED) {
39 42 : for (int i = 0; i < features.size(); i++) {
40 21 : features[i]->onDetection(candidate);
41 : }
42 :
43 21 : if (detectionAction.valid()) {
44 6 : if (clone)
45 0 : detectionAction->process(candidate->clone(false));
46 : else
47 6 : detectionAction->process(candidate);
48 : }
49 :
50 21 : if (!flagKey.empty())
51 5 : candidate->setProperty(flagKey, flagValue);
52 :
53 21 : if (makeInactive)
54 18 : candidate->setActive(false);
55 : }
56 476 : }
57 :
58 2 : void Observer::setFlag(std::string key, std::string value) {
59 2 : flagKey = key;
60 2 : flagValue = value;
61 2 : }
62 :
63 0 : std::string Observer::getDescription() const {
64 0 : std::stringstream ss;
65 0 : ss << "Observer";
66 0 : for (int i = 0; i < features.size(); i++)
67 0 : ss << "\n " << features[i]->getDescription() << "\n";
68 0 : ss << " Flag: '" << flagKey << "' -> '" << flagValue << "'\n";
69 0 : ss << " MakeInactive: " << (makeInactive ? "yes\n" : "no\n");
70 0 : if (detectionAction.valid())
71 0 : ss << " Action: " << detectionAction->getDescription() << ", clone: " << (clone ? "yes" : "no");
72 :
73 0 : return ss.str();
74 0 : }
75 :
76 6 : void Observer::setDeactivateOnDetection(bool deactivate) {
77 6 : makeInactive = deactivate;
78 6 : }
79 :
80 : // ObserverFeature ------------------------------------------------------------
81 0 : DetectionState ObserverFeature::checkDetection(Candidate *candidate) const {
82 0 : return NOTHING;
83 : }
84 :
85 21 : void ObserverFeature::onDetection(Candidate *candidate) const {
86 21 : }
87 :
88 0 : std::string ObserverFeature::getDescription() const {
89 0 : return description;
90 : }
91 :
92 : // ObserverDetectAll ----------------------------------------------------------
93 2 : DetectionState ObserverDetectAll::checkDetection(Candidate *candidate) const {
94 2 : return DETECTED;
95 : }
96 :
97 0 : std::string ObserverDetectAll::getDescription() const {
98 0 : return description;
99 : }
100 :
101 : // ObserverTracking --------------------------------------------------------
102 0 : ObserverTracking::ObserverTracking(Vector3d center, double radius, double stepSize) :
103 0 : center(center), radius(radius), stepSize(stepSize) {
104 : if (stepSize == 0) {
105 : stepSize = radius / 10.;
106 : }
107 0 : }
108 :
109 0 : DetectionState ObserverTracking::checkDetection(Candidate *candidate) const {
110 : // current distance to observer sphere center
111 0 : double d = (candidate->current.getPosition() - center).getR();
112 :
113 : // no detection if outside of observer sphere
114 0 : if (d > radius) {
115 : // conservatively limit next step to prevent overshooting
116 0 : candidate->limitNextStep(fabs(d - radius));
117 :
118 0 : return NOTHING;
119 : } else {
120 : // limit next step
121 0 : candidate->limitNextStep(stepSize);
122 :
123 0 : return DETECTED;
124 : }
125 : }
126 :
127 0 : std::string ObserverTracking::getDescription() const {
128 0 : std::stringstream ss;
129 0 : ss << "ObserverTracking: ";
130 0 : ss << "center = " << center / Mpc << " Mpc, ";
131 0 : ss << "radius = " << radius / Mpc << " Mpc";
132 0 : ss << "stepSize = " << stepSize / Mpc << " Mpc";
133 0 : return ss.str();
134 0 : }
135 :
136 : // Observer1D --------------------------------------------------------------
137 456 : DetectionState Observer1D::checkDetection(Candidate *candidate) const {
138 456 : double x = candidate->current.getPosition().x;
139 456 : if (x > 0) {
140 : // Limits the next step size to prevent candidates from overshooting in case of non-detection
141 449 : candidate->limitNextStep(x);
142 449 : return NOTHING;
143 : }
144 : // Detects particles when reaching x = 0
145 : return DETECTED;
146 : }
147 :
148 0 : std::string Observer1D::getDescription() const {
149 0 : return "Observer1D: observer at x = 0";
150 : }
151 :
152 : // ObserverRedshiftWindow -----------------------------------------------------
153 0 : ObserverRedshiftWindow::ObserverRedshiftWindow(double zmin, double zmax) :
154 0 : zmin(zmin), zmax(zmax) {
155 0 : }
156 :
157 0 : DetectionState ObserverRedshiftWindow::checkDetection(
158 : Candidate *candidate) const {
159 0 : double z = candidate->getRedshift();
160 0 : if (z > zmax)
161 : return VETO;
162 0 : if (z < zmin)
163 0 : return VETO;
164 : return NOTHING;
165 : }
166 :
167 0 : std::string ObserverRedshiftWindow::getDescription() const {
168 0 : std::stringstream ss;
169 0 : ss << "ObserverRedshiftWindow: z = " << zmin << " - " << zmax;
170 0 : return ss.str();
171 0 : }
172 :
173 : // ObserverInactiveVeto -------------------------------------------------------
174 0 : DetectionState ObserverInactiveVeto::checkDetection(Candidate *c) const {
175 0 : if (not(c->isActive()))
176 0 : return VETO;
177 : return NOTHING;
178 : }
179 :
180 0 : std::string ObserverInactiveVeto::getDescription() const {
181 0 : return "ObserverInactiveVeto";
182 : }
183 :
184 : // ObserverNucleusVeto --------------------------------------------------------
185 0 : DetectionState ObserverNucleusVeto::checkDetection(Candidate *c) const {
186 0 : if (isNucleus(c->current.getId()))
187 0 : return VETO;
188 : return NOTHING;
189 : }
190 :
191 0 : std::string ObserverNucleusVeto::getDescription() const {
192 0 : return "ObserverNucleusVeto";
193 : }
194 :
195 : // ObserverNeutrinoVeto -------------------------------------------------------
196 0 : DetectionState ObserverNeutrinoVeto::checkDetection(Candidate *c) const {
197 0 : int id = abs(c->current.getId());
198 0 : if ((id == 12) or (id == 14) or (id == 16))
199 0 : return VETO;
200 : return NOTHING;
201 : }
202 :
203 0 : std::string ObserverNeutrinoVeto::getDescription() const {
204 0 : return "ObserverNeutrinoVeto";
205 : }
206 :
207 : // ObserverPhotonVeto ---------------------------------------------------------
208 0 : DetectionState ObserverPhotonVeto::checkDetection(Candidate *c) const {
209 0 : if (c->current.getId() == 22)
210 0 : return VETO;
211 : return NOTHING;
212 : }
213 :
214 0 : std::string ObserverPhotonVeto::getDescription() const {
215 0 : return "ObserverPhotonVeto";
216 : }
217 :
218 : // ObserverElectronVeto ---------------------------------------------------------
219 0 : DetectionState ObserverElectronVeto::checkDetection(Candidate *c) const {
220 0 : if (abs(c->current.getId()) == 11)
221 0 : return VETO;
222 : return NOTHING;
223 : }
224 :
225 0 : std::string ObserverElectronVeto::getDescription() const {
226 0 : return "ObserverElectronVeto";
227 : }
228 :
229 : // ObserverCustomVeto -------------------------------------------------------
230 0 : ObserverParticleIdVeto::ObserverParticleIdVeto(int pId) {
231 0 : vetoParticleId = pId;
232 0 : }
233 :
234 0 : DetectionState ObserverParticleIdVeto::checkDetection(Candidate *c) const {
235 0 : if (c->current.getId() == vetoParticleId)
236 0 : return VETO;
237 : return NOTHING;
238 : }
239 :
240 0 : std::string ObserverParticleIdVeto::getDescription() const {
241 0 : return "ObserverParticleIdVeto";
242 : }
243 :
244 :
245 : // ObserverTimeEvolution --------------------------------------------------------
246 0 : ObserverTimeEvolution::ObserverTimeEvolution() {}
247 :
248 1 : ObserverTimeEvolution::ObserverTimeEvolution(double min, double dist, double numb) {
249 : setIsLogarithmicScaling(false);
250 1 : setMaximum(min + (numb - 1) * dist);
251 1 : setMinimum(min);
252 1 : setNIntervals(numb);
253 1 : }
254 :
255 2 : ObserverTimeEvolution::ObserverTimeEvolution(double min, double max, double numb, bool log) {
256 : setIsLogarithmicScaling(log);
257 2 : setMinimum(min);
258 : setMaximum(max);
259 2 : setNIntervals(numb);
260 2 : }
261 :
262 1 : ObserverTimeEvolution::ObserverTimeEvolution(const std::vector<double> &detList){
263 1 : setTimes(detList);
264 1 : }
265 :
266 9 : DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
267 :
268 9 : if (nIntervals) {
269 9 : double length = c->getTrajectoryLength();
270 : size_t index;
271 9 : const std::string DI = "DetectionIndex";
272 :
273 : // Load the last detection index
274 9 : if (c->hasProperty(DI)) {
275 10 : index = c->getProperty(DI).asUInt64();
276 : }
277 : else {
278 : index = 0;
279 : }
280 :
281 : // Break if the particle has been detected once for all possible times.
282 9 : if (index >= nIntervals) {
283 : return NOTHING;
284 : }
285 :
286 : // Calculate the distance to next detection
287 9 : double distance = length - getTime(index);
288 :
289 : // Limit next step and detect candidate.
290 : // Increase the index by one in case of detection
291 9 : if (distance < 0.) {
292 4 : c->limitNextStep(-distance);
293 : return NOTHING;
294 : }
295 : else {
296 :
297 5 : if (index < nIntervals-2) {
298 1 : c->limitNextStep(getTime(index+1)-length);
299 : }
300 5 : c->setProperty(DI, Variant::fromUInt64(index+1));
301 :
302 5 : return DETECTED;
303 : }
304 : }
305 : return NOTHING;
306 : }
307 :
308 3 : void ObserverTimeEvolution::clear(){
309 3 : doDetListConstruction = false;
310 : detList.clear();
311 3 : detList.resize(0);
312 : setNIntervals(0);
313 3 : }
314 :
315 12 : void ObserverTimeEvolution::constructDetListIfEmpty(){
316 12 : if (detList.empty() && doDetListConstruction) {
317 : std::vector<double> detListTemp;
318 : size_t counter = 0;
319 7 : while (getTime(counter)<=maximum) {
320 6 : detListTemp.push_back(getTime(counter));
321 6 : counter++;
322 : }
323 1 : detList.assign(detListTemp.begin(), detListTemp.end());
324 1 : }
325 12 : }
326 :
327 12 : void ObserverTimeEvolution::addTime(const double &time){
328 12 : constructDetListIfEmpty();
329 12 : detList.push_back(time);
330 12 : setNIntervals(nIntervals + 1); // increase number of entries by one
331 12 : }
332 :
333 2 : void ObserverTimeEvolution::addTimeRange(double min, double max, double numb, bool log) {
334 12 : for (size_t i = 0; i < numb; i++) {
335 10 : if (log) {
336 4 : if ( min <= 0 ){
337 : std::cout << "min can not be <= 0 if log=true" << std::endl;
338 0 : throw new std::runtime_error("min can not be <= 0 if log=true");
339 : }
340 4 : addTime(min * pow(max / min, i / (numb - 1.0)));
341 : } else {
342 6 : addTime(min + i * (max - min) / (numb - 1.0));
343 : }
344 : }
345 : // allready corrected by addTime, just here to be safe
346 2 : setNIntervals(detList.size());
347 2 : }
348 :
349 1 : void ObserverTimeEvolution::setTimes(const std::vector<double> &detList){
350 1 : this->detList.assign(detList.begin(), detList.end());
351 1 : setNIntervals(detList.size());
352 1 : setMinimum(detList.front());
353 1 : setMaximum(detList.back());
354 1 : doDetListConstruction = false;
355 1 : }
356 :
357 4 : void ObserverTimeEvolution::setMinimum(double min){
358 4 : if ( (min <= 0) && isLogarithmicScaling){
359 : std::cout << "minimum can not be <= 0 if isLogarithmicScaling=true" << std::endl;
360 0 : throw new std::runtime_error("minimum can not be <= 0 if isLogarithmicScaling=true");
361 : }
362 4 : this->minimum = min;
363 4 : }
364 :
365 65 : double ObserverTimeEvolution::getTime(size_t index) const {
366 65 : if (!detList.empty()) {
367 36 : return detList.at(index);
368 29 : } else if (isLogarithmicScaling) {
369 6 : return minimum * pow(maximum / minimum, index / (nIntervals - 1.0));
370 : } else {
371 23 : return minimum + index * (maximum - minimum) / (nIntervals - 1.0);
372 : }
373 : }
374 :
375 9 : const std::vector<double>& ObserverTimeEvolution::getTimes() const {
376 9 : tempDetList.resize(nIntervals);
377 51 : for (size_t i = 0; i < nIntervals; i++){
378 42 : tempDetList[i] = getTime(i);
379 : }
380 9 : return tempDetList;
381 : }
382 :
383 0 : std::string ObserverTimeEvolution::getDescription() const {
384 0 : std::stringstream s;
385 0 : s << "List of Detection lengths in kpc";
386 0 : for (size_t i = 0; i < nIntervals; i++)
387 0 : s << " - " << getTime(i) / kpc;
388 0 : return s.str();
389 0 : }
390 :
391 : // ObserverSurface--------------------------------------------------------------
392 2 : ObserverSurface::ObserverSurface(Surface* _surface) : surface(_surface) { }
393 :
394 4 : DetectionState ObserverSurface::checkDetection(Candidate *candidate) const
395 : {
396 4 : double currentDistance = surface->distance(candidate->current.getPosition());
397 4 : double previousDistance = surface->distance(candidate->previous.getPosition());
398 4 : candidate->limitNextStep(fabs(currentDistance));
399 :
400 4 : if (currentDistance * previousDistance > 0)
401 : return NOTHING;
402 2 : else if (previousDistance == 0)
403 : return NOTHING;
404 : else
405 2 : return DETECTED;
406 : }
407 :
408 0 : std::string ObserverSurface::getDescription() const {
409 0 : std::stringstream ss;
410 0 : ss << "ObserverSurface: << " << surface->getDescription();
411 0 : return ss.str();
412 0 : }
413 :
414 :
415 : } // namespace crpropa
|
