Line data    Source code

       1             : #include "crpropa/module/BreakCondition.h"
       2             : #include "crpropa/ParticleID.h"
       3             : #include "crpropa/Units.h"
       4             : 
       5             : #include <sstream>
       6             : 
       7             : namespace crpropa {
       8             : 
       9          13 : MaximumTrajectoryLength::MaximumTrajectoryLength(double maxLength) :
      10          13 :                 maxLength(maxLength) {
      11          13 : }
      12             : 
      13           0 : void MaximumTrajectoryLength::setMaximumTrajectoryLength(double length) {
      14           0 :         maxLength = length;
      15           0 : }
      16             : 
      17           0 : double MaximumTrajectoryLength::getMaximumTrajectoryLength() const {
      18           0 :         return maxLength;
      19             : }
      20             : 
      21           1 : void MaximumTrajectoryLength::addObserverPosition(const Vector3d& position) {
      22           1 :         observerPositions.push_back(position);
      23           1 : }
      24             : 
      25           0 : const std::vector<Vector3d>& MaximumTrajectoryLength::getObserverPositions() const {
      26           0 :         return observerPositions;
      27             : }
      28             : 
      29           0 : std::string MaximumTrajectoryLength::getDescription() const {
      30           0 :         std::stringstream s;
      31           0 :         s << "Maximum trajectory length: " << maxLength / Mpc << " Mpc, ";
      32           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
      33           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
      34           0 :         if (rejectAction.valid())
      35           0 :                 s << ", Action: " << rejectAction->getDescription();
      36           0 :         s << "\n  Observer positions: \n";
      37           0 :         for (size_t i = 0; i < observerPositions.size(); i++)
      38           0 :                 s << "    - " << observerPositions[i] / Mpc << " Mpc\n";
      39           0 :         return s.str();
      40           0 : }
      41             : 
      42      491649 : void MaximumTrajectoryLength::process(Candidate *c) const {
      43      491649 :         double length = c->getTrajectoryLength();
      44      491649 :         Vector3d position = c->current.getPosition();
      45             : 
      46      491649 :         if(observerPositions.size()) {
      47             :                 bool inRange = false;
      48           4 :                 for (size_t i = 0; i < observerPositions.size(); i++) {
      49           2 :                         double distance = position.getDistanceTo(observerPositions[i]);
      50           2 :                         if (distance + length < maxLength)
      51             :                                 inRange = true;
      52             :                 }
      53           2 :                 if (!inRange) {
      54           1 :                         reject(c);
      55             :                         return;
      56             :                 }
      57             :         }
      58             : 
      59      491648 :         if (length >= maxLength) {
      60        1107 :                 reject(c);
      61             :         } else {
      62      490541 :                 c->limitNextStep(maxLength - length);
      63             :         }
      64             : }
      65             : 
      66             : //*****************************************************************************
      67           2 : MinimumEnergy::MinimumEnergy(double minEnergy) :
      68           2 :                 minEnergy(minEnergy) {
      69           2 : }
      70             : 
      71           0 : void MinimumEnergy::setMinimumEnergy(double energy) {
      72           0 :         minEnergy = energy;
      73           0 : }
      74             : 
      75           0 : double MinimumEnergy::getMinimumEnergy() const {
      76           0 :         return minEnergy;
      77             : }
      78             : 
      79       32602 : void MinimumEnergy::process(Candidate *c) const {
      80       32602 :         if (c->current.getEnergy() > minEnergy)
      81             :                 return;
      82             :         else
      83           1 :                 reject(c);
      84             : }
      85             : 
      86           0 : std::string MinimumEnergy::getDescription() const {
      87           0 :         std::stringstream s;
      88           0 :         s << "Minimum energy: " << minEnergy / EeV << " EeV, ";
      89           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
      90           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
      91           0 :         if (rejectAction.valid())
      92           0 :                 s << ", Action: " << rejectAction->getDescription();
      93           0 :         return s.str();
      94           0 : }
      95             : 
      96             : //*****************************************************************************
      97           0 : MinimumRigidity::MinimumRigidity(double minRigidity) :
      98           0 :                 minRigidity(minRigidity) {
      99           0 : }
     100             : 
     101           0 : void MinimumRigidity::setMinimumRigidity(double minRigidity) {
     102           0 :         this->minRigidity = minRigidity;
     103           0 : }
     104             : 
     105           0 : double MinimumRigidity::getMinimumRigidity() const {
     106           0 :         return minRigidity;
     107             : }
     108             : 
     109           0 : void MinimumRigidity::process(Candidate *c) const {
     110           0 :         if (c->current.getRigidity() < minRigidity)
     111           0 :                 reject(c);
     112           0 : }
     113             : 
     114           0 : std::string MinimumRigidity::getDescription() const {
     115           0 :         std::stringstream s;
     116           0 :         s << "Minimum rigidity: " << minRigidity / EeV << " EeV, ";
     117           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     118           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     119           0 :         if (rejectAction.valid())
     120           0 :                 s << ", Action: " << rejectAction->getDescription();
     121           0 :         return s.str();
     122           0 : }
     123             : 
     124             : //*****************************************************************************
     125           1 : MinimumRedshift::MinimumRedshift(double zmin) :
     126           1 :                 zmin(zmin) {
     127           1 : }
     128             : 
     129           0 : void MinimumRedshift::setMinimumRedshift(double z) {
     130           0 :         zmin = z;
     131           0 : }
     132             : 
     133           0 : double MinimumRedshift::getMinimumRedshift() {
     134           0 :         return zmin;
     135             : }
     136             : 
     137           2 : void MinimumRedshift::process(Candidate* c) const {
     138           2 :         if (c->getRedshift() > zmin)
     139             :                 return;
     140             :         else
     141           1 :                 reject(c);
     142             : }
     143             : 
     144           0 : std::string MinimumRedshift::getDescription() const {
     145           0 :         std::stringstream s;
     146           0 :         s << "Minimum redshift: " << zmin << ", ";
     147           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     148           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     149           0 :         if (rejectAction.valid())
     150           0 :                 s << ", Action: " << rejectAction->getDescription();
     151           0 :         return s.str();
     152           0 : }
     153             : 
     154             : //*****************************************************************************
     155           1 : MinimumChargeNumber::MinimumChargeNumber(int minChargeNumber) :
     156           1 :                 minChargeNumber(minChargeNumber) {
     157           1 : }
     158             : 
     159           0 : void MinimumChargeNumber::setMinimumChargeNumber(int chargeNumber) {
     160           0 :         minChargeNumber = chargeNumber;
     161           0 : }
     162             : 
     163           0 : int MinimumChargeNumber::getMinimumChargeNumber() const {
     164           0 :         return minChargeNumber;
     165             : }
     166             : 
     167           4 : void MinimumChargeNumber::process(Candidate *c) const {
     168           4 :         if (chargeNumber(c->current.getId()) > minChargeNumber)
     169             :                 return;
     170             :         else
     171           2 :                 reject(c);
     172             : }
     173             : 
     174           0 : std::string MinimumChargeNumber::getDescription() const {
     175           0 :         std::stringstream s;
     176           0 :         s << "Minimum charge number: " << minChargeNumber;
     177           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     178           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     179           0 :         if (rejectAction.valid())
     180           0 :                 s << ", Action: " << rejectAction->getDescription();
     181           0 :         return s.str();
     182           0 : }
     183             : 
     184             : //*****************************************************************************
     185           1 : MinimumEnergyPerParticleId::MinimumEnergyPerParticleId(double minEnergyOthers) {
     186           1 :         setMinimumEnergyOthers(minEnergyOthers);
     187           1 : }
     188             : 
     189           3 : void MinimumEnergyPerParticleId::add(int id, double energy) {
     190           3 :         particleIds.push_back(id);
     191           3 :         minEnergies.push_back(energy);
     192           3 : }
     193             : 
     194           1 : void MinimumEnergyPerParticleId::setMinimumEnergyOthers(double energy) {
     195           1 :         minEnergyOthers = energy;
     196           1 : }
     197             : 
     198           0 : double MinimumEnergyPerParticleId::getMinimumEnergyOthers() const {
     199           0 :         return minEnergyOthers;
     200             : }
     201             : 
     202           5 : void MinimumEnergyPerParticleId::process(Candidate *c) const {
     203          15 :         for (int i = 0; i < particleIds.size(); i++) {
     204          12 :                 if (c->current.getId() == particleIds[i]) {
     205           5 :                         if (c->current.getEnergy() < minEnergies[i])
     206           3 :                                 reject(c);
     207             :                         else
     208             :                                 return;
     209             :                 }
     210             :         }
     211             : 
     212           3 :         if (c->current.getEnergy() < minEnergyOthers)
     213           1 :                 reject(c);
     214             :         else
     215             :                 return;
     216             : }
     217             : 
     218           0 : std::string MinimumEnergyPerParticleId::getDescription() const {
     219           0 :         std::stringstream s;
     220           0 :         s << "Minimum energy for non-specified particles: " << minEnergyOthers / eV << " eV";
     221           0 :         for (int i = 0; i < minEnergies.size(); i++) {
     222           0 :                 s << "  for particle " << particleIds[i] << " : " << minEnergies[i] / eV << " eV";
     223             :         }
     224           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     225           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     226           0 :         if (rejectAction.valid())
     227           0 :                 s << ", Action: " << rejectAction->getDescription();
     228           0 :         return s.str();
     229           0 : }
     230             : 
     231             : //*****************************************************************************
     232           1 : DetectionLength::DetectionLength(double detLength) :
     233           1 :                 detLength(detLength) {
     234           1 : }
     235             : 
     236           0 : void DetectionLength::setDetectionLength(double length) {
     237           0 :         detLength = length;
     238           0 : }
     239             : 
     240           0 : double DetectionLength::getDetectionLength() const {
     241           0 :         return detLength;
     242             : }
     243             : 
     244             : 
     245           0 : std::string DetectionLength::getDescription() const {
     246           0 :         std::stringstream s;
     247           0 :         s << "Detection length: " << detLength / kpc << " kpc, ";
     248           0 :         s << "Flag: '" << rejectFlagKey << "' -> '" << rejectFlagValue << "', ";
     249           0 :         s << "MakeInactive: " << (makeRejectedInactive ? "yes" : "no");
     250           0 :         if (rejectAction.valid())
     251           0 :                 s << ", Action: " << rejectAction->getDescription();
     252           0 :         return s.str();
     253           0 : }
     254             : 
     255           2 : void DetectionLength::process(Candidate *c) const {
     256           2 :         double length = c->getTrajectoryLength();
     257           2 :         double step = c->getCurrentStep();
     258             : 
     259           2 :         if (length >= detLength && length - step < detLength) {
     260           1 :                 reject(c);
     261             :         } else {
     262           1 :                 c->limitNextStep(detLength - length);
     263             :         }
     264           2 : }
     265             : 
     266             : 
     267             : } // namespace crpropa