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       32639 : void Observer::process(Candidate *candidate) const {
      28             :         // loop over all features and have them check the particle
      29             :         DetectionState state = NOTHING;
      30       65278 :         for (int i = 0; i < features.size(); i++) {
      31       32639 :                 DetectionState s = features[i]->checkDetection(candidate);
      32       32639 :                 if (s == VETO)
      33             :                         state = VETO;
      34       32639 :                 else if ((s == DETECTED) && (state != VETO))
      35             :                         state = DETECTED;
      36             :         }
      37             : 
      38       32639 :         if (state == DETECTED) {
      39         146 :                 for (int i = 0; i < features.size(); i++) {
      40          73 :                         features[i]->onDetection(candidate);
      41             :                 }
      42             : 
      43          73 :                 if (detectionAction.valid()) {
      44          59 :                         if (clone)
      45           0 :                                 detectionAction->process(candidate->clone(false));
      46             :                         else
      47          59 :                                 detectionAction->process(candidate);
      48             :                 }
      49             : 
      50          73 :                 if (!flagKey.empty())
      51           4 :                         candidate->setProperty(flagKey, flagValue);
      52             : 
      53          73 :                 if (makeInactive)
      54          69 :                         candidate->setActive(false);
      55             :         }
      56       32639 : }
      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           2 : void Observer::setDeactivateOnDetection(bool deactivate) {
      77           2 :         makeInactive = deactivate;
      78           2 : }
      79             : 
      80             : // ObserverFeature ------------------------------------------------------------
      81           0 : DetectionState ObserverFeature::checkDetection(Candidate *candidate) const {
      82           0 :         return NOTHING;
      83             : }
      84             : 
      85          73 : void ObserverFeature::onDetection(Candidate *candidate) const {
      86          73 : }
      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       32622 : DetectionState Observer1D::checkDetection(Candidate *candidate) const {
     138       32622 :         double x = candidate->current.getPosition().x;
     139       32622 :         if (x > 0) {
     140             :                 // Limits the next step size to prevent candidates from overshooting in case of non-detection
     141       32562 :                 candidate->limitNextStep(x);
     142       32562 :                 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           1 :         double max = min + numb * dist;
     250             :         bool log = false;
     251           1 :         addTimeRange(min, max, numb, log);
     252           1 : }
     253             : 
     254           1 : ObserverTimeEvolution::ObserverTimeEvolution(double min, double max, double numb, bool log) {
     255           1 :         addTimeRange(min, max, numb, log);
     256           1 : }
     257             : 
     258             : 
     259           6 : DetectionState ObserverTimeEvolution::checkDetection(Candidate *c) const {
     260             : 
     261           6 :         if (detList.size()) {
     262           6 :                 double length = c->getTrajectoryLength();
     263             :                 size_t index;
     264           6 :                 const std::string DI = "DetectionIndex";
     265             :                 std::string value;
     266             : 
     267             :                 // Load the last detection index
     268           6 :                 if (c->hasProperty(DI)) {
     269           4 :                         index = c->getProperty(DI).asUInt64();
     270             :                 }
     271             :                 else {
     272             :                         index = 0;
     273             :                 }
     274             : 
     275             :                 // Break if the particle has been detected once for all detList entries.
     276           6 :                 if (index > detList.size()) {
     277             :                         return NOTHING;
     278             :                 }
     279             : 
     280             :                 // Calculate the distance to next detection
     281           6 :                 double distance = length - detList[index];
     282             : 
     283             :                 // Limit next step and detect candidate.
     284             :                 // Increase the index by one in case of detection
     285           6 :                 if (distance < 0.) {
     286           2 :                         c->limitNextStep(-distance);
     287             :                         return NOTHING;
     288             :                 }
     289             :                 else {
     290             : 
     291           4 :                         if (index < detList.size()-1) {
     292           2 :                                 c->limitNextStep(detList[index+1]-length);
     293             :                         }
     294           4 :                         c->setProperty(DI, Variant::fromUInt64(index+1));
     295             : 
     296           4 :                         return DETECTED;
     297             :                 }
     298             :         }
     299             :         return NOTHING;
     300             : }
     301             : 
     302           4 : void ObserverTimeEvolution::addTime(const double& t) {
     303           4 :         detList.push_back(t);
     304           4 : }
     305             : 
     306           2 : void ObserverTimeEvolution::addTimeRange(double min, double max, double numb, bool log) {
     307           6 :         for (size_t i = 0; i < numb; i++) {
     308           4 :                 if (log == true) {
     309           2 :                         addTime(min * pow(max / min, i / (numb - 1.0)));
     310             :                 } else {
     311           2 :                         addTime(min + i * (max - min) / numb);
     312             :                 }
     313             :         }
     314           2 : }
     315             : 
     316           0 : const std::vector<double>& ObserverTimeEvolution::getTimes() const {
     317           0 :         return detList;
     318             : }
     319             : 
     320           0 : std::string ObserverTimeEvolution::getDescription() const {
     321           0 :         std::stringstream s;
     322           0 :         s << "List of Detection lengths in kpc";
     323           0 :         for (size_t i = 0; i < detList.size(); i++)
     324           0 :           s << "  - " << detList[i] / kpc;
     325           0 :         return s.str();
     326           0 : }
     327             : 
     328             : // ObserverSurface--------------------------------------------------------------
     329           2 : ObserverSurface::ObserverSurface(Surface* _surface) : surface(_surface) { }
     330             : 
     331           4 : DetectionState ObserverSurface::checkDetection(Candidate *candidate) const
     332             : {
     333           4 :                 double currentDistance = surface->distance(candidate->current.getPosition());
     334           4 :                 double previousDistance = surface->distance(candidate->previous.getPosition());
     335           4 :                 candidate->limitNextStep(fabs(currentDistance));
     336             : 
     337           4 :                 if (currentDistance * previousDistance > 0)
     338             :                         return NOTHING;
     339           2 :                 else if (previousDistance == 0)
     340             :                         return NOTHING;
     341             :                 else
     342           2 :                         return DETECTED;
     343             : }
     344             : 
     345           0 : std::string ObserverSurface::getDescription() const {
     346           0 :         std::stringstream ss;
     347           0 :         ss << "ObserverSurface: << " << surface->getDescription();
     348           0 :         return ss.str();
     349           0 : }
     350             : 
     351             : } // namespace crpropa