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