Line data    Source code

       1             : #include "crpropa/Candidate.h"
       2             : #include "crpropa/ParticleID.h"
       3             : #include "crpropa/Units.h"
       4             : 
       5             : #include <stdexcept>
       6             : 
       7             : namespace crpropa {
       8             : 
       9     3371882 : Candidate::Candidate(int id, double E, Vector3d pos, Vector3d dir, double z, double weight, std::string tagOrigin) :
      10     3371882 :   redshift(z), trajectoryLength(0), weight(weight), currentStep(0), nextStep(0), active(true), parent(0), tagOrigin(tagOrigin) {
      11     3371882 :         ParticleState state(id, E, pos, dir);
      12             :         source = state;
      13             :         created = state;
      14             :         previous = state;
      15             :         current = state;
      16             : 
      17             : #if defined(OPENMP_3_1)
      18             :                 #pragma omp atomic capture
      19             :                 {serialNumber = nextSerialNumber++;}
      20             : #elif defined(__GNUC__)
      21     3371882 :                 {serialNumber = __sync_add_and_fetch(&nextSerialNumber, 1);}
      22             : #else
      23             :                 #pragma omp critical
      24             :                 {serialNumber = nextSerialNumber++;}
      25             : #endif
      26             : 
      27     3371882 : }
      28             : 
      29          20 : Candidate::Candidate(const ParticleState &state) :
      30          20 :                 source(state), created(state), current(state), previous(state), redshift(0), trajectoryLength(0), currentStep(0), nextStep(0), active(true), parent(0), tagOrigin ("PRIM") {
      31             : 
      32             : #if defined(OPENMP_3_1)
      33             :                 #pragma omp atomic capture
      34             :                 {serialNumber = nextSerialNumber++;}
      35             : #elif defined(__GNUC__)
      36          20 :                 {serialNumber = __sync_add_and_fetch(&nextSerialNumber, 1);}
      37             : #else
      38             :                 #pragma omp critical
      39             :                 {serialNumber = nextSerialNumber++;}
      40             : #endif
      41             : 
      42          20 : }
      43             : 
      44       46315 : bool Candidate::isActive() const {
      45       46315 :         return active;
      46             : }
      47             : 
      48        1759 : void Candidate::setActive(bool b) {
      49        1759 :         active = b;
      50        1759 : }
      51             : 
      52      679017 : double Candidate::getRedshift() const {
      53      679017 :         return redshift;
      54             : }
      55             : 
      56     1011739 : double Candidate::getTrajectoryLength() const {
      57     1011739 :         return trajectoryLength;
      58             : }
      59             : 
      60        1085 : double Candidate::getWeight() const {
      61        1085 :         return weight;
      62             : }
      63             : 
      64      262738 : double Candidate::getCurrentStep() const {
      65      262738 :         return currentStep;
      66             : }
      67             : 
      68      524343 : double Candidate::getNextStep() const {
      69      524343 :         return nextStep;
      70             : }
      71             : 
      72     3360758 : void Candidate::setRedshift(double z) {
      73     3360758 :         redshift = z;
      74     3360758 : }
      75             : 
      76     3327970 : void Candidate::setTrajectoryLength(double a) {
      77     3327970 :         trajectoryLength = a;
      78     3327970 : }
      79             : 
      80     3328954 : void Candidate::setWeight(double w) {
      81     3328954 :         weight = w;
      82     3328954 : }
      83             : 
      84           5 : void Candidate::updateWeight(double w) {
      85           5 :   weight *= w;
      86           5 : }
      87             : 
      88      526034 : void Candidate::setCurrentStep(double lstep) {
      89      526034 :         currentStep = lstep;
      90      526034 :         trajectoryLength += lstep;
      91      526034 : }
      92             : 
      93      524358 : void Candidate::setNextStep(double step) {
      94      524358 :         nextStep = step;
      95      524358 : }
      96             : 
      97      654514 : void Candidate::limitNextStep(double step) {
      98      654514 :         nextStep = std::min(nextStep, step);
      99      654514 : }
     100             : 
     101        1143 : void Candidate::setProperty(const std::string &name, const Variant &value) {
     102        1143 :         properties[name] = value;
     103        1143 : }
     104             : 
     105     3327943 : void Candidate::setTagOrigin (std::string tagOrigin) {
     106     3327943 :         this->tagOrigin = tagOrigin;
     107     3327943 : }
     108             : 
     109          82 : std::string Candidate::getTagOrigin () const {
     110          82 :         return tagOrigin;
     111             : }
     112             : 
     113          11 : const Variant &Candidate::getProperty(const std::string &name) const {
     114          11 :         PropertyMap::const_iterator i = properties.find(name);
     115          11 :         if (i == properties.end())
     116           0 :                 throw std::runtime_error("Unknown candidate property: " + name);
     117          11 :         return i->second;
     118             : }
     119             : 
     120           5 : bool Candidate::removeProperty(const std::string& name) {
     121           5 :         PropertyMap::iterator i = properties.find(name);
     122           5 :         if (i == properties.end())
     123             :                 return false;
     124             :         properties.erase(i);
     125             :         return true;
     126             : }
     127             : 
     128          30 : bool Candidate::hasProperty(const std::string &name) const {
     129          30 :         PropertyMap::const_iterator i = properties.find(name);
     130          30 :         if (i == properties.end())
     131          10 :                 return false;
     132             :         return true;
     133             : }
     134             : 
     135           4 : void Candidate::addSecondary(Candidate *c) {
     136           4 :         secondaries.push_back(c);
     137           4 : }
     138             : 
     139           2 : void Candidate::addSecondary(int id, double energy, double w, std::string tagOrigin) {
     140           4 :         ref_ptr<Candidate> secondary = new Candidate;
     141           2 :         secondary->setRedshift(redshift);
     142           2 :         secondary->setTrajectoryLength(trajectoryLength);
     143           2 :         secondary->setWeight(weight * w);
     144           4 :         secondary->setTagOrigin(tagOrigin);
     145           2 :         for (PropertyMap::const_iterator it = properties.begin(); it != properties.end(); ++it) {
     146           0 :                 secondary->setProperty(it->first, it->second);         
     147             :         }
     148             :         secondary->source = source;
     149             :         secondary->previous = previous;
     150             :         secondary->created = previous;
     151           2 :         secondary->current = current;
     152           2 :         secondary->current.setId(id);
     153           2 :         secondary->current.setEnergy(energy);
     154           2 :         secondary->parent = this;
     155           2 :         secondaries.push_back(secondary);
     156           2 : }
     157             : 
     158     3327939 : void Candidate::addSecondary(int id, double energy, Vector3d position, double w, std::string tagOrigin) {
     159     6655878 :         ref_ptr<Candidate> secondary = new Candidate;
     160     3327939 :         secondary->setRedshift(redshift);
     161     3327939 :         secondary->setTrajectoryLength(trajectoryLength - (current.getPosition() - position).getR());
     162     3327939 :         secondary->setWeight(weight * w);
     163     6655878 :         secondary->setTagOrigin(tagOrigin);
     164     3327939 :         for (PropertyMap::const_iterator it = properties.begin(); it != properties.end(); ++it) {
     165           0 :                 secondary->setProperty(it->first, it->second);         
     166             :         }
     167             :         secondary->source = source;
     168             :         secondary->previous = previous;
     169     3327939 :         secondary->created = previous;
     170     3327939 :         secondary->current = current;
     171     3327939 :         secondary->current.setId(id);
     172     3327939 :         secondary->current.setEnergy(energy);
     173     3327939 :         secondary->current.setPosition(position);
     174     3327939 :         secondary->created.setPosition(position);
     175     3327939 :         secondary->parent = this;
     176     3327939 :         secondaries.push_back(secondary);
     177     3327939 : }
     178             : 
     179           0 : void Candidate::clearSecondaries() {
     180           0 :         secondaries.clear();
     181           0 : }
     182             : 
     183           0 : std::string Candidate::getDescription() const {
     184           0 :         std::stringstream ss;
     185           0 :         ss << "CosmicRay at z = " << getRedshift() << "\n";
     186           0 :         ss << "  source:  " << source.getDescription() << "\n";
     187           0 :         ss << "  current: " << current.getDescription();
     188           0 :         return ss.str();
     189           0 : }
     190             : 
     191          30 : ref_ptr<Candidate> Candidate::clone(bool recursive) const {
     192          60 :         ref_ptr<Candidate> cloned = new Candidate;
     193             :         cloned->source = source;
     194             :         cloned->created = created;
     195             :         cloned->current = current;
     196             :         cloned->previous = previous;
     197             : 
     198          30 :         cloned->properties = properties;
     199          30 :         cloned->active = active;
     200          30 :         cloned->redshift = redshift;
     201          30 :         cloned->weight = weight;
     202          30 :         cloned->trajectoryLength = trajectoryLength;
     203          30 :         cloned->currentStep = currentStep;
     204          30 :         cloned->nextStep = nextStep;
     205          30 :         if (recursive) {
     206           0 :                 cloned->secondaries.reserve(secondaries.size());
     207           0 :                 for (size_t i = 0; i < secondaries.size(); i++) {
     208           0 :                         ref_ptr<Candidate> s = secondaries[i]->clone(recursive);
     209           0 :                         s->parent = cloned;
     210           0 :                         cloned->secondaries.push_back(s);
     211             :                 }
     212             :         }
     213          30 :         return cloned;
     214             : }
     215             : 
     216         126 : uint64_t Candidate::getSerialNumber() const {
     217         126 :         return serialNumber;
     218             : }
     219             : 
     220          11 : void Candidate::setSerialNumber(const uint64_t snr) {
     221          11 :         serialNumber = snr;
     222          11 : }
     223             : 
     224          70 : uint64_t Candidate::getSourceSerialNumber() const {
     225         142 :         if (parent)
     226             :                 return parent->getSourceSerialNumber();
     227             :         else
     228          70 :                 return serialNumber;
     229             : }
     230             : 
     231          69 : uint64_t Candidate::getCreatedSerialNumber() const {
     232          69 :         if (parent)
     233          56 :                 return parent->getSerialNumber();
     234             :         else
     235          13 :                 return serialNumber;
     236             : }
     237             : 
     238           1 : void Candidate::setNextSerialNumber(uint64_t snr) {
     239           1 :         nextSerialNumber = snr;
     240           1 : }
     241             : 
     242           6 : uint64_t Candidate::getNextSerialNumber() {
     243           6 :         return nextSerialNumber;
     244             : }
     245             : 
     246             : uint64_t Candidate::nextSerialNumber = 0;
     247             : 
     248           1 : void Candidate::restart() {
     249           1 :         setActive(true);
     250           1 :         setTrajectoryLength(0);
     251             :         previous = source;
     252             :         current = source;
     253           1 : }
     254             : 
     255             : } // namespace crpropa