Line data Source code
1 : #ifndef CRPROPA_SOURCE_H
2 : #define CRPROPA_SOURCE_H
3 :
4 : #include "crpropa/Candidate.h"
5 : #include "crpropa/Grid.h"
6 : #include "crpropa/EmissionMap.h"
7 : #include "crpropa/massDistribution/Density.h"
8 :
9 :
10 : #include <vector>
11 :
12 : namespace crpropa {
13 : /** @addtogroup SourceFeatures
14 : * @{
15 : */
16 :
17 :
18 : /**
19 : @class SourceFeature
20 : @brief Abstract base class for specific source features
21 : */
22 7 : class SourceFeature: public Referenced {
23 : protected:
24 : std::string description;
25 : public:
26 0 : virtual void prepareParticle(ParticleState& particle) const {};
27 : virtual void prepareCandidate(Candidate& candidate) const;
28 : std::string getDescription() const;
29 : };
30 :
31 :
32 : /**
33 : @class SourceInterface
34 : @brief Abstract base class for sources
35 : */
36 : class SourceInterface : public Referenced {
37 : public:
38 : virtual ref_ptr<Candidate> getCandidate() const = 0;
39 : virtual std::string getDescription() const = 0;
40 : };
41 :
42 :
43 : /**
44 : @class Source
45 : @brief General source of particles
46 :
47 : This class is a container for source features.
48 : The source prepares a new candidate by passing it to all its source features
49 : to be modified accordingly.
50 : */
51 10 : class Source: public SourceInterface {
52 : std::vector<ref_ptr<SourceFeature> > features;
53 : public:
54 : void add(SourceFeature* feature);
55 : ref_ptr<Candidate> getCandidate() const;
56 : std::string getDescription() const;
57 : };
58 :
59 :
60 : /**
61 : @class SourceList
62 : @brief List of particle sources of individual luminosities.
63 :
64 : The SourceList is a source itself. It can be used if several sources are
65 : needed in one simulation.
66 : */
67 3 : class SourceList: public SourceInterface {
68 : std::vector<ref_ptr<Source> > sources;
69 : std::vector<double> cdf;
70 : public:
71 : /** Add an individual source to the list.
72 : @param source source to be added
73 : @param weight weight of the source; defaults to 1.
74 : */
75 : void add(Source* source, double weight = 1);
76 : ref_ptr<Candidate> getCandidate() const;
77 : std::string getDescription() const;
78 : };
79 :
80 :
81 : /**
82 : @class SourceParticleType
83 : @brief Particle type at the source
84 :
85 : This feature assigns a single particle type to the source.
86 : For multiple types, use e.g. SourceMultipleParticleTypes.
87 : Particles are identified following the PDG numbering scheme:
88 : https://pdg.lbl.gov/2019/reviews/rpp2019-rev-monte-carlo-numbering.pdf
89 : */
90 : class SourceParticleType: public SourceFeature {
91 : int id;
92 : public:
93 : /** Constructor for a source with a sign
94 : @param id id of the particle following the PDG numbering scheme
95 : */
96 : SourceParticleType(int id);
97 : void prepareParticle(ParticleState &particle) const;
98 : void setDescription();
99 : };
100 :
101 :
102 : /**
103 : @class SourceMultipleParticleTypes
104 : @brief Multiple particle types with individual relative abundances
105 :
106 : This feature assigns particle types to the events emitted by the sources.
107 : It is possible to control the relative abundance of each particle species.
108 : Particles are identified following the PDG numbering scheme:
109 : https://pdg.lbl.gov/2019/reviews/rpp2019-rev-monte-carlo-numbering.pdf
110 : */
111 : class SourceMultipleParticleTypes: public SourceFeature {
112 : std::vector<int> particleTypes;
113 : std::vector<double> cdf;
114 : public:
115 : /** Constructor
116 : */
117 : SourceMultipleParticleTypes();
118 : /** Add an individual particle type.
119 : @param id id of the particle following the PDG numbering scheme
120 : @param weight relative abundance of individual particle species
121 : */
122 : void add(int id, double weight = 1);
123 : void prepareParticle(ParticleState &particle) const;
124 : void setDescription();
125 : };
126 :
127 :
128 : /**
129 : @class SourceEnergy
130 : @brief Sets the initial energy of the emitted particles to a specific value
131 :
132 : This feature assigns a monochromatic spectrum, i.e., a single energy to all particles.
133 : */
134 : class SourceEnergy: public SourceFeature {
135 : double E;
136 : public:
137 : /** Constructor
138 : @param energy energy of the particle (in Joules)
139 : */
140 : SourceEnergy(double energy);
141 : void prepareParticle(ParticleState &particle) const;
142 : void setDescription();
143 : };
144 :
145 :
146 : /**
147 : @class SourcePowerLawSpectrum
148 : @brief Particle energy following a power-law spectrum
149 :
150 : The power law is of the form: dN/dE ~ E^index, for energies in the interval [Emin, Emax].
151 : */
152 1 : class SourcePowerLawSpectrum: public SourceFeature {
153 : double Emin;
154 : double Emax;
155 : double index;
156 : public:
157 : /** Constructor
158 : @param Emin minimum energy (in Joules)
159 : @param Emax maximum energy (in Joules)
160 : @param index spectral index of the power law
161 : */
162 : SourcePowerLawSpectrum(double Emin, double Emax, double index);
163 : void prepareParticle(ParticleState &particle) const;
164 : void setDescription();
165 : };
166 :
167 :
168 : /**
169 : @class SourceComposition
170 : @brief Multiple nuclear species with a rigidity-dependent power-law spectrum
171 :
172 : The power law is of the form: E^index, for energies in the interval [Emin, Z * Rmax].
173 : */
174 : class SourceComposition: public SourceFeature {
175 : double Emin;
176 : double Rmax;
177 : double index;
178 : std::vector<int> nuclei;
179 : std::vector<double> cdf;
180 : public:
181 : /** Constructor
182 : @param Emin minimum energy (in Joules)
183 : @param Rmax maximum rigidity (in Volts)
184 : @param index spectral index of the power law
185 : */
186 : SourceComposition(double Emin, double Rmax, double index);
187 : /** Add individual particle species with a given abundance
188 : @param id id of the particle following the PDG numbering scheme
189 : @param abundance relative abundance of the particle species
190 : */
191 : void add(int id, double abundance);
192 : /** Add individual particle species with a given abundance
193 : @param A atomic mass of the cosmic-ray nucleus
194 : @param Z atomic number of the cosmic-ray nucleus
195 : @param abundance relative abundance of the particle species
196 : */
197 : void add(int A, int Z, double abundance);
198 : void prepareParticle(ParticleState &particle) const;
199 : void setDescription();
200 : };
201 :
202 :
203 : /**
204 : @class SourcePosition
205 : @brief Position of a point source
206 : */
207 1 : class SourcePosition: public SourceFeature {
208 : Vector3d position;
209 : public:
210 : /** Constructor for a source in 3D
211 : @param position vector containing the coordinates of the point source [in meters]
212 : */
213 : SourcePosition(Vector3d position);
214 : /** Constructor for a source in 1D
215 : @param d distance of the point source to the observer at x = 0 [in meters];
216 : internally this will be converted to a vector with x-coordinate equal to d
217 : */
218 : SourcePosition(double d);
219 : void prepareParticle(ParticleState &state) const;
220 : void setDescription();
221 : };
222 :
223 :
224 : /**
225 : @class SourceMultiplePositions
226 : @brief Multiple point-source positions with individual luminosities
227 : */
228 : class SourceMultiplePositions: public SourceFeature {
229 : std::vector<Vector3d> positions;
230 : std::vector<double> cdf;
231 : public:
232 : /** Constructor.
233 : The sources must be added individually to the object.
234 : */
235 : SourceMultiplePositions();
236 : /** Add an individual source with a given luminosity/contribution.
237 : @param position vector containing the coordinates of the point source [in meters]
238 : @param weight luminosity/contribution of the individual source
239 : */
240 : void add(Vector3d position, double weight = 1);
241 : void prepareParticle(ParticleState &particle) const;
242 : void setDescription();
243 : };
244 :
245 :
246 : /**
247 : @class SourceUniformSphere
248 : @brief Uniform distribution of sources in a spherical volume
249 : */
250 1 : class SourceUniformSphere: public SourceFeature {
251 : Vector3d center;
252 : double radius;
253 : public:
254 : /** Constructor
255 : @param center vector containing the coordinates of the center of the sphere
256 : @param radius radius of the sphere
257 : */
258 : SourceUniformSphere(Vector3d center, double radius);
259 : void prepareParticle(ParticleState &particle) const;
260 : void setDescription();
261 : };
262 :
263 :
264 : /**
265 : @class SourceUniformHollowSphere
266 : @brief Uniform distribution of sources between two spheres
267 : */
268 1 : class SourceUniformHollowSphere: public SourceFeature {
269 : Vector3d center;
270 : double radius_inner;
271 : double radius_outer;
272 : public:
273 : /** Constructor
274 : @param center vector containing the coordinates of the center of the sphere
275 : @param radius_inner radius of the inner sphere
276 : @param radius_outer radius of the outer sphere
277 : */
278 : SourceUniformHollowSphere(Vector3d center,
279 : double radius_inner, double radius_outer);
280 : void prepareParticle(ParticleState &particle) const;
281 : void setDescription();
282 : };
283 :
284 :
285 : /**
286 : @class SourceUniformShell
287 : @brief Uniform distribution of source positions on the surface of a sphere
288 : */
289 : class SourceUniformShell: public SourceFeature {
290 : Vector3d center;
291 : double radius;
292 : public:
293 : /** Constructor
294 : @param center vector containing the coordinates of the center of the sphere
295 : @param radius radius of the sphere
296 : */
297 : SourceUniformShell(Vector3d center, double radius);
298 : void prepareParticle(ParticleState &particle) const;
299 : void setDescription();
300 : };
301 :
302 :
303 : /**
304 : @class SourceUniformBox
305 : @brief Uniform random source positions inside a box. The box is aligned with the coordinate axes.
306 : */
307 1 : class SourceUniformBox: public SourceFeature {
308 : Vector3d origin; // lower box corner
309 : Vector3d size; // sizes along each coordinate axes.
310 : public:
311 : /** Constructor
312 : @param origin vector corresponding to the lower box corner
313 : @param size vector corresponding to the box sizes along each direction
314 : */
315 : SourceUniformBox(Vector3d origin, Vector3d size);
316 : void prepareParticle(ParticleState &particle) const;
317 : void setDescription();
318 : };
319 :
320 :
321 : /**
322 : @class SourceUniformCylinder
323 : @brief Uniform distribution of source positions inside the volume of a cylinder whose axis is along the z-axis.
324 :
325 : The circle of the cylinder lays in the xy-plane and the height is along the z-axis.
326 : */
327 1 : class SourceUniformCylinder: public SourceFeature {
328 : Vector3d origin; // central point of cylinder
329 : double height; // total height of the cylinder along z-axis. Half over/under the center.
330 : double radius; // radius of the cylinder in the xy-plane
331 : public:
332 : /** Constructor
333 : @param origin vector corresponding to the center of the cylinder axis
334 : @param height height of the cylinder, half lays over the origin, half is lower
335 : @param radius radius of the cylinder
336 : */
337 : SourceUniformCylinder(Vector3d origin, double height, double radius);
338 : void prepareParticle(ParticleState &particle) const;
339 : void setDescription();
340 : };
341 :
342 :
343 : /**
344 : @class SourceSNRDistribution
345 : @brief Source distribution that follows the Galactic SNR distribution in 2D
346 :
347 : The origin of the distribution is the Galactic center. The default maximum radius is set
348 : to rMax=20 kpc and the default maximum height is zMax = 5 kpc.
349 : See G. Case and D. Bhattacharya (1996) for the details of the distribution.
350 : */
351 1 : class SourceSNRDistribution: public SourceFeature {
352 : double rEarth; // parameter given by observation
353 : double alpha; // parameter to shift the maximum in R direction
354 : double beta; // parameter to shift the maximum in R direction
355 : double zg; // exponential cut parameter in z direction
356 : double frMax; // helper for efficient sampling
357 : double fzMax; // helper for efficient sampling
358 : double rMax; // maximum radial distance - default 20 kpc
359 : // (due to the extension of the JF12 field)
360 : double zMax; // maximum distance from galactic plane - default 5 kpc
361 : void setFrMax(); // calculate frMax with the current parameter.
362 :
363 : public:
364 : /** Default constructor.
365 : Default parameters are:
366 : . rEarth = 8.5 kpc
367 : . alpha = 2
368 : . beta = 3.53
369 : . zg = 300 pc
370 : . rMax = 20 kpc
371 : . zMax = 5 kpc
372 : */
373 : SourceSNRDistribution();
374 : /** Generic constructor
375 : @param rEarth distance from Earth to the Galactic centre [in meters]
376 : @param alpha parameter that shifts radially the maximum of the distributions
377 : @param beta parameter that shifts radially the maximum of the distributions
378 : @param zg exponential cut-off parameter in the z-direction [in meters]
379 : */
380 : SourceSNRDistribution(double rEarth,double alpha, double beta, double zg);
381 :
382 : void prepareParticle(ParticleState &particle) const;
383 : /**
384 : radial distribution of the SNR density.
385 : @param r galactocentric radius in [meter]
386 : */
387 : double fr(double r) const;
388 : /**
389 : height distribution of the SNR density.
390 : @param z height over/under the galactic plane in [meter]
391 : */
392 : double fz(double z) const;
393 :
394 : /**
395 : Set the exponential cut-off parameter in the z-direction.
396 : @param Zg cut-off parameter
397 : */
398 : void setFzMax(double Zg);
399 :
400 : /**
401 : @param rMax maximal radius up to which sources are possible
402 : */
403 : void setRMax(double rMax);
404 :
405 : /**
406 : @param zMax maximal height up to which sources are possible
407 : */
408 : void setZMax(double zMax);
409 :
410 : // parameter for the raidal distribution
411 : void setAlpha(double a);
412 : // parameter for the exponential cut-off in the radial distribution
413 : void setBeta(double b);
414 : double getFrMax() const;
415 : double getFzMax() const;
416 : double getRMax() const;
417 : double getZMax() const;
418 : double getAlpha() const;
419 : double getBeta() const;
420 : void setDescription();
421 : };
422 :
423 :
424 : /**
425 : @class SourcePulsarDistribution
426 : @brief Source distribution following the Galactic pulsar distribution
427 :
428 : A logarithmic spiral with four arms is used for the radial distribution.
429 : The z-distribution is a simple exponentially decaying distribution.
430 : The pulsar distribution is explained in detail in C.-A. Faucher-Giguere
431 : and V. M. Kaspi, ApJ 643 (May, 2006) 332. The radial distribution is
432 : parametrized as in Blasi and Amato, JCAP 1 (Jan., 2012) 10.
433 : */
434 : class SourcePulsarDistribution: public SourceFeature {
435 : double rEarth; // parameter given by observation
436 : double beta; // parameter to shift the maximum in R direction
437 : double zg; // exponential cut parameter in z direction
438 : double frMax; // helper for efficient sampling
439 : double fzMax; // helper for efficient sampling
440 : double rMax; // maximum radial distance - default 22 kpc
441 : double zMax; // maximum distance from galactic plane - default 5 kpc
442 : double rBlur; // relative smearing factor for the radius
443 : double thetaBlur; // smearing factor for the angle. Unit = [1/length]
444 : public:
445 : /** Default constructor.
446 : Default parameters are:
447 : . rEarth = 8.5 kpc
448 : . beta = 3.53
449 : . zg = 300 pc
450 : . Rmax = 22 kpc
451 : . Zmax = 5 kpc
452 : . rBlur = 0.07
453 : . thetaBlur = 0.35 / kpc
454 : */
455 : SourcePulsarDistribution();
456 : /** Generic constructor
457 : @param rEarth distance from Earth to the Galactic centre [in meters]
458 : @param beta parameter that shifts radially the maximum of the distributions
459 : @param zg exponential cut-off parameter in the z-direction [in meters]
460 : @param rBlur relative smearing factor for radius
461 : @param thetaBlur smearing factor for the angle [in 1 / meters]
462 : */
463 : SourcePulsarDistribution(double rEarth, double beta, double zg, double rBlur, double thetaBlur);
464 : void prepareParticle(ParticleState &particle) const;
465 :
466 : /**
467 : radial distribution of pulsars
468 : @param r galactocentric radius
469 : */
470 : double fr(double r) const;
471 : /**
472 : z distribution of pulsars
473 : @param z height over/under the galactic plane
474 : */
475 : double fz(double z) const;
476 : double ftheta(int i, double r) const;
477 : double blurR(double r_tilde) const;
478 : double blurTheta(double theta_tilde, double r_tilde) const;
479 : void setFrMax(double R, double b);
480 : void setFzMax(double zg);
481 : void setRMax(double rMax);
482 : void setZMax(double zMax);
483 : void setRBlur(double rBlur);
484 : void setThetaBlur(double thetaBlur);
485 : double getFrMax();
486 : double getFzMax();
487 : double getRMax();
488 : double getZMax();
489 : double getRBlur();
490 : double getThetaBlur();
491 : void setDescription();
492 : };
493 :
494 :
495 : /**
496 : @class SourceUniform1D
497 : @brief Uniform source distribution in 1D
498 :
499 : This source property sets random x-coordinates according to a uniform source
500 : distribution in a given distance interval. If cosmological effects are included,
501 : this is done by drawing a light-travel distance from a flat distribution and
502 : converting to a comoving distance. In the absence of cosmological effects, the
503 : positions are drawn uniformly in the light-travel distance interval (as opposed
504 : to a comoving interval).
505 : The source positions are assigned to the x-coordinate (Vector3d(distance, 0, 0))
506 : in this one-dimensional case.
507 : */
508 : class SourceUniform1D: public SourceFeature {
509 : double minD; // minimum light-travel distance
510 : double maxD; // maximum light-travel distance
511 : bool withCosmology; // whether to account for cosmological effects (expansion of the Universe)
512 : public:
513 : /** Constructor
514 : @param minD minimum distance; comoving if withCosmology is True
515 : @param maxD maximum distance; comoving if withCosmology is True
516 : @param withCosmology whether to account for cosmological effects (expansion of the Universe)
517 : */
518 : SourceUniform1D(double minD, double maxD, bool withCosmology = true);
519 : void prepareParticle(ParticleState& particle) const;
520 : void setDescription();
521 : };
522 :
523 :
524 : /**
525 : @class SourceDensityGrid
526 : @brief Random source positions from a density grid
527 : */
528 : class SourceDensityGrid: public SourceFeature {
529 : ref_ptr<Grid1f> grid;
530 : public:
531 : /** Constructor
532 : @param densityGrid 3D grid containing the density of sources in each cell
533 : */
534 : SourceDensityGrid(ref_ptr<Grid1f> densityGrid);
535 : void prepareParticle(ParticleState &particle) const;
536 : void setDescription();
537 : };
538 :
539 :
540 : /**
541 : @class SourceDensityGrid1D
542 : @brief Random source positions from a 1D density grid
543 : */
544 : class SourceDensityGrid1D: public SourceFeature {
545 : ref_ptr<Grid1f> grid; // 1D grid with Ny = Nz = 1
546 : public:
547 : /** Constructor
548 : @param densityGrid 1D grid containing the density of sources in each cell, Ny and Nz must be 1
549 : */
550 : SourceDensityGrid1D(ref_ptr<Grid1f> densityGrid);
551 : void prepareParticle(ParticleState &particle) const;
552 : void setDescription();
553 : };
554 :
555 :
556 : /**
557 : @class SourceIsotropicEmission
558 : @brief Isotropic emission from a source
559 : */
560 : class SourceIsotropicEmission: public SourceFeature {
561 : public:
562 : /** Constructor
563 : */
564 : SourceIsotropicEmission();
565 : void prepareParticle(ParticleState &particle) const;
566 : void setDescription();
567 : };
568 :
569 :
570 : /**
571 : @class SourceDirectedEmission
572 : @brief Directed emission from a source from the von-Mises-Fisher distribution
573 :
574 : The emission from the source is generated following the von-Mises-Fisher distribution
575 : with mean direction mu and concentration parameter kappa.
576 : The sampling from the vMF distribution follows this document by Julian Straub:
577 : http://people.csail.mit.edu/jstraub/download/straub2017vonMisesFisherInference.pdf
578 : The emitted particles are assigned a weight so that the detected particles can be
579 : reweighted to an isotropic emission distribution instead of a vMF distribution.
580 : For details, see PoS (ICRC2019) 447.
581 : */
582 1 : class SourceDirectedEmission: public SourceFeature {
583 : Vector3d mu; // Mean emission direction in the vMF distribution
584 : double kappa; // Concentration parameter of the vMF distribution
585 : public:
586 : /** Constructor
587 : @param mu mean direction of the emission, mu should be normelized
588 : @param kappa concentration parameter
589 : */
590 : SourceDirectedEmission(Vector3d mu, double kappa);
591 : void prepareCandidate(Candidate &candidate) const;
592 : void setDescription();
593 : };
594 :
595 : /**
596 : @class SourceLambertDistributionOnSphere
597 : @brief Uniform random position on a sphere with isotropic Lamberts distributed directions.
598 :
599 : This function should be used for crosschecking the arrival distribution for a
600 : Galactic propagation with an isotropic arrival distribution at the Edge of our
601 : Galaxy. Note, that for simulation speed you should rather use the backtracking
602 : technique: see e.g. http://physik.rwth-aachen.de/parsec
603 : */
604 : class SourceLambertDistributionOnSphere: public SourceFeature {
605 : Vector3d center; // center of the sphere
606 : double radius; // radius of the sphere
607 : bool inward; // if true, direction point inwards
608 : public:
609 : /** Constructor
610 : @param center vector containing the coordinates of the center of the sphere
611 : @param radius radius of the sphere
612 : @param inward if true, the directions point inwards
613 : */
614 : SourceLambertDistributionOnSphere(const Vector3d ¢er, double radius, bool inward);
615 : void prepareParticle(ParticleState &particle) const;
616 : void setDescription();
617 : };
618 :
619 :
620 : /**
621 : @class SourceDirection
622 : @brief Collimated emission along a specific direction
623 : */
624 : class SourceDirection: public SourceFeature {
625 : Vector3d direction;
626 : public:
627 : /** Constructor
628 : @param direction Vector3d corresponding to the direction of emission
629 : */
630 : SourceDirection(Vector3d direction = Vector3d(-1, 0, 0));
631 : void prepareParticle(ParticleState &particle) const;
632 : void setDescription();
633 : };
634 :
635 :
636 : /**
637 : @class SourceEmissionMap
638 : @brief Deactivate Candidate if it has zero probability in provided EmissionMap.
639 :
640 : This feature does not change the direction of the candidate. Therefore a usefull direction feature (isotropic or directed emission)
641 : must be added to the sources before. The propability of the emission map is not taken into account.
642 : */
643 : class SourceEmissionMap: public SourceFeature {
644 : ref_ptr<EmissionMap> emissionMap;
645 : public:
646 : /** Constructor
647 : @param emissionMap emission map containing probabilities of emission in various directions
648 : */
649 : SourceEmissionMap(EmissionMap *emissionMap);
650 : void prepareCandidate(Candidate &candidate) const;
651 : void setEmissionMap(EmissionMap *emissionMap);
652 : void setDescription();
653 : };
654 :
655 :
656 : /**
657 : @class SourceEmissionCone
658 : @brief Uniform emission within a cone
659 : */
660 1 : class SourceEmissionCone: public SourceFeature {
661 : Vector3d direction;
662 : double aperture;
663 : public:
664 : /** Constructor
665 : @param direction Vector3d corresponding to the cone axis
666 : @param aperture opening angle of the cone
667 : */
668 : SourceEmissionCone(Vector3d direction, double aperture);
669 : void prepareParticle(ParticleState &particle) const;
670 :
671 : /**
672 : @param direction Vector3d corresponding to the cone axis
673 : */
674 : void setDirection(Vector3d direction);
675 : void setDescription();
676 : };
677 :
678 :
679 : /**
680 : @class SourceRedshift
681 : @brief Emission of particles at a specific redshift (or time)
682 :
683 : The redshift coordinate is used to treat cosmological effects and as a time coordinate.
684 : Consider, for instance, a source located at a distance corresponding to a redshift z.
685 : In the absence of processes that cause time delays (e.g., magnetic deflections), particles
686 : from this source could arrive after a time corresponding to the source redshift. Charged
687 : particles, on the other hand, can arrive at a time later than the corresponding straight-
688 : line travel duration.
689 : This treatment is also useful for time-dependent studies (e.g. transient sources).
690 : */
691 : class SourceRedshift: public SourceFeature {
692 : double z;
693 : public:
694 : /** Constructor
695 : @param z redshift of emission
696 : */
697 : SourceRedshift(double z);
698 : void prepareCandidate(Candidate &candidate) const;
699 : void setDescription();
700 : };
701 :
702 :
703 : /**
704 : @class SourceUniformRedshift
705 : @brief Random redshift (time of emission) from a uniform distribution
706 :
707 : This function assigns random redshifts to the particles emitted by a given source.
708 : These values are drawn from a uniform redshift distribution in the interval [zmin, zmax].
709 : The redshift coordinate is used to treat cosmological effects and as a time coordinate.
710 : Consider, for instance, a source located at a distance corresponding to a redshift z.
711 : In the absence of processes that cause time delays (e.g., magnetic deflections), particles
712 : from this source could arrive after a time corresponding to the source redshift. Charged
713 : particles, on the other hand, can arrive at a time later than the corresponding straight-
714 : line travel duration.
715 : This treatment is also useful for time-dependent studies (e.g. transient sources).
716 : */
717 : class SourceUniformRedshift: public SourceFeature {
718 : double zmin, zmax;
719 : public:
720 : /** Constructor
721 : @param zmin minimum redshift
722 : @param zmax maximum redshift
723 : */
724 : SourceUniformRedshift(double zmin, double zmax);
725 : void prepareCandidate(Candidate &candidate) const;
726 : void setDescription();
727 : };
728 :
729 :
730 : /**
731 : @class SourceRedshiftEvolution
732 : @brief Random redshift (time of emission) from (1+z)^m distribution
733 :
734 : This assigns redshifts to a given source according to a typical power-law distribution.
735 : The redshift coordinate is used to treat cosmological effects and as a time coordinate.
736 : Consider, for instance, a source located at a distance corresponding to a redshift z.
737 : In the absence of processes that cause time delays (e.g., magnetic deflections), particles
738 : from this source could arrive after a time corresponding to the source redshift. Charged
739 : particles, on the other hand, can arrive at a time later than the corresponding straight-
740 : line travel duration.
741 : This treatment is also useful for time-dependent studies (e.g. transient sources).
742 : */
743 2 : class SourceRedshiftEvolution: public SourceFeature {
744 : double zmin, zmax;
745 : double m;
746 : public:
747 : /** Constructor
748 : @param m index of the power law (1 + z)^m
749 : @param zmin minimum redshift
750 : @param zmax maximum redshift
751 : */
752 : SourceRedshiftEvolution(double m, double zmin, double zmax);
753 : void prepareCandidate(Candidate &candidate) const;
754 : };
755 :
756 :
757 : /**
758 : @class SourceRedshift1D
759 : @brief Redshift according to the distance to 0
760 :
761 : This source property sets the redshift according to the distance from
762 : the source to the origin (0, 0, 0).
763 : It must be added after the position of the source is set because it
764 : computes the redshifts based on the source distance.
765 : */
766 : class SourceRedshift1D: public SourceFeature {
767 : public:
768 : /** Constructor
769 : */
770 : SourceRedshift1D();
771 : void prepareCandidate(Candidate &candidate) const;
772 : void setDescription();
773 : };
774 :
775 :
776 : #ifdef CRPROPA_HAVE_MUPARSER
777 : /**
778 : @class SourceGenericComposition
779 : @brief Add multiple cosmic rays with energies described by an expression string
780 :
781 : This is particularly useful if an arbitrary combination of nuclei types with
782 : specific energy spectra. The strings parsed may contain 'A' (atomic mass),
783 : 'Z' (atomic number). The following units are recognized as part of the strings:
784 : GeV, TeV, PeV, EeV. The variable for energy is 'E', with limits 'Emin', 'Emax'.
785 : This property only works if muparser is available.
786 : For details about the library see:
787 : https://beltoforion.de/en/muparser/
788 : */
789 : class SourceGenericComposition: public SourceFeature {
790 : public:
791 7 : struct Nucleus {
792 : int id;
793 : std::vector<double> cdf;
794 : };
795 : /** Constructor
796 : @param Emin minimum energy [in Joules]
797 : @param Emax maximum energy [in Joules]
798 : @param expression string containing the expression to generate the composition
799 : @param bins number of energy bins
800 : */
801 : SourceGenericComposition(double Emin, double Emax, std::string expression, size_t bins = 1024);
802 : /** Add an individual particle id.
803 : @param id id of the particle following the PDG numbering scheme
804 : @param abundance relative abundance of individual particle species
805 : */
806 : void add(int id, double abundance);
807 : /** Add an individual particle id.
808 : @param A atomic mass of the cosmic-ray nucleus
809 : @param Z atomic number of the cosmic-ray nucleus
810 : @param abundance relative abundance of individual particle species
811 : */
812 : void add(int A, int Z, double abundance);
813 : void prepareParticle(ParticleState &particle) const;
814 : void setDescription();
815 :
816 : const std::vector<double> *getNucleusCDF(int id) const {
817 0 : for (size_t i = 0; i < nuclei.size(); i++) {
818 0 : if (nuclei[i].id == id)
819 0 : return &nuclei[i].cdf;
820 : }
821 : return 0;
822 : }
823 :
824 : protected:
825 : double Emin, Emax;
826 : size_t bins;
827 : std::string expression;
828 : std::vector<double> energy;
829 :
830 : std::vector<Nucleus> nuclei;
831 : std::vector<double> cdf;
832 :
833 : };
834 : #endif
835 :
836 : /**
837 : * @class SourceTag
838 : * @brief All candidates from this source get a given tag. This can be used to distinguish between different sources that follow the same spatial distribution
839 : *
840 : * Sets the tag of the candidate. Can be used to trace back additional candidate properties, e.g. production interaction or source type.
841 : * The interaction overwrites the candidate tag from the source for all secondaries.
842 : */
843 :
844 : class SourceTag: public SourceFeature {
845 : private:
846 : std::string sourceTag;
847 :
848 : public:
849 : SourceTag(std::string tag);
850 : void prepareCandidate(Candidate &candidate) const;
851 : void setDescription();
852 : void setTag(std::string tag);
853 : };
854 :
855 : /**
856 : @class SourceMassDistribution
857 : @brief Source position follows a given mass distribution
858 :
859 : The (source)position of the candidate is sampled from a given mass distribution. The distribution uses the getDensity function of the density module.
860 : If a weighting for different components is desired, the use of different densities in a densityList is recommended.
861 :
862 : The sampling range of the position can be restricted. Default is a sampling for x in [-20, 20] * kpc, y in [-20, 20] * kpc and z in [-4, 4] * kpc.
863 : */
864 : class SourceMassDistribution: public SourceFeature {
865 : private:
866 : ref_ptr<Density> density; //< density distribution
867 : double maxDensity; //< maximal value of the density in the region of interest
868 : double xMin, xMax; //< x-range to sample positions
869 : double yMin, yMax; //< y-range to sample positions
870 : double zMin, zMax; //< z-range to sample positions
871 : int maxTries = 10000; //< maximal number of tries to sample the position
872 :
873 : public:
874 : /** Constructor
875 : @param density: CRPropa mass distribution
876 : @param maxDensity: maximal density in the region where the position should be sampled
877 : @param x: the position will be sampled in the range [-x, x]. Non symmetric values can be set with setXrange.
878 : @param y: the position will be sampled in the range [-y, y]. Non symmetric values can be set with setYrange.
879 : @param z: the position will be sampled in the range [-z, z]. Non symmetric values can be set with setZrange.
880 : */
881 : SourceMassDistribution(ref_ptr<Density> density, double maxDensity = 0, double x = 20 * kpc, double y = 20 * kpc, double z = 4 * kpc);
882 :
883 : void prepareParticle(ParticleState &particle) const;
884 :
885 : /** Set the maximal density in the region of interest. This parameter is necessary for the sampling
886 : @param maxDensity: maximal density in [particle / m^3]
887 : */
888 : void setMaximalDensity(double maxDensity);
889 :
890 : /** set x-range in which the position of the candidate will be sampled. x in [xMin, xMax].
891 : @param xMin: minimal x value of the allowed sample range in [m]
892 : @param xMax: maximal x value of the allowed sample range in [m]
893 : */
894 : void setXrange(double xMin, double xMax);
895 :
896 : /** set y-range in which the position of the candidate will be sampled. y in [yMin, yMax].
897 : @param yMin: minimal y value of the allowed sample range in [m]
898 : @param yMax: maximal y value of the allowed sample range in [m]
899 : */
900 : void setYrange(double yMin, double yMax);
901 :
902 : /** set z-range in which the position of the candidate will be sampled. z in [zMin, zMax].
903 : @param zMin: minimal z value of the allowed sample range in [m]
904 : @param zMax: maximal z value of the allowed sample range in [m]
905 : */
906 : void setZrange(double zMin, double zMax);
907 :
908 : /* samples the position. Can be used for testing.
909 : @return Vector3d with sampled position
910 : */
911 : Vector3d samplePosition() const;
912 :
913 :
914 : /** set the number of maximal tries until the sampling routine breaks.
915 : @param tries: number of the maximal tries
916 : */
917 : void setMaximalTries(int tries);
918 :
919 : std::string getDescription();
920 : };
921 :
922 : /** @} */ // end of group SourceFeature
923 :
924 : } // namespace crpropa
925 :
926 : #endif // CRPROPA_SOURCE_H
|
