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 436 : void MinimumEnergy::process(Candidate *c) const {
80 436 : 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
|
