Line data Source code
1 : #include "crpropa/module/EMTripletPairProduction.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/Random.h"
4 : #include "crpropa/Common.h"
5 :
6 : #include <vector>
7 : #include <cmath>
8 :
9 : namespace crpropa {
10 :
11 : static const double mec2 = mass_electron * c_squared;
12 :
13 3 : EMTripletPairProduction::EMTripletPairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit, ref_ptr<Surface> surface) {
14 6 : setSurface(surface);
15 3 : setPhotonField(photonField);
16 3 : setHaveElectrons(haveElectrons);
17 3 : setLimit(limit);
18 3 : setThinning(thinning);
19 3 : }
20 :
21 15 : void EMTripletPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
22 15 : this->photonField = photonField;
23 15 : std::string fname = photonField->getFieldName();
24 30 : setDescription("EMTripletPairProduction: " + fname);
25 15 : if (!this->photonField->hasPositionDependence()){
26 : this->interactionRates = new InteractionRatesHomogeneous(
27 45 : getDataPath("EMTripletPairProduction/rate_" + fname + ".txt"),
28 45 : getDataPath("EMTripletPairProduction/cdf_" + fname + ".txt")
29 30 : );
30 : } else {
31 : this->interactionRates = new InteractionRatesPositionDependent(
32 0 : getDataPath("EMTripletPairProduction/"+fname+"/Rate/"),
33 0 : getDataPath("EMTripletPairProduction/"+fname+"/CumulativeRate/"),
34 : this->surface
35 0 : );
36 : }
37 15 : }
38 :
39 4 : void EMTripletPairProduction::setHaveElectrons(bool haveElectrons) {
40 4 : this->haveElectrons = haveElectrons;
41 4 : }
42 :
43 3 : void EMTripletPairProduction::setLimit(double limit) {
44 3 : this->limit = limit;
45 3 : }
46 :
47 3 : void EMTripletPairProduction::setThinning(double thinning) {
48 3 : this->thinning = thinning;
49 3 : }
50 :
51 3 : void EMTripletPairProduction::setSurface(ref_ptr<Surface> surface) {
52 3 : this->surface = surface;
53 3 : }
54 :
55 0 : ref_ptr<Surface> EMTripletPairProduction::getSurface() const {
56 0 : return this->surface;
57 : }
58 :
59 0 : void EMTripletPairProduction::setInteractionRates(ref_ptr<InteractionRates> intRates) {
60 0 : this->interactionRates = intRates;
61 0 : }
62 :
63 0 : ref_ptr<InteractionRates> EMTripletPairProduction::getInteractionRates() const {
64 0 : return this->interactionRates;
65 : }
66 :
67 0 : void EMTripletPairProduction::initRate(std::string path) {
68 0 : this->interactionRates->initRate(path);
69 0 : }
70 :
71 0 : void EMTripletPairProduction::initCumulativeRate(std::string path) {
72 0 : this->interactionRates->initCumulativeRate(path);
73 0 : }
74 :
75 1 : double EMTripletPairProduction::getRate(double E, const Vector3d &position, double z) const {
76 1 : return this->interactionRates->getProcessRate(E, position) * pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
77 : }
78 :
79 1 : void EMTripletPairProduction::performInteraction(Candidate *candidate) const {
80 :
81 : // scale the particle energy instead of background photons
82 1 : double z = candidate->getRedshift();
83 1 : double E = candidate->current.getEnergy() * (1 + z);
84 1 : Vector3d position = candidate->current.getPosition();
85 :
86 : std::vector<double> tabE;
87 : std::vector<double> tabs;
88 : std::vector<std::vector<double>> tabCDF;
89 :
90 1 : this->interactionRates->loadPerformInteractionTabs(position, tabE, tabs, tabCDF);
91 :
92 1 : if (E < tabE.front() or E > tabE.back())
93 : return;
94 :
95 : // sample the value of eps
96 1 : Random &random = Random::instance();
97 1 : size_t i = closestIndex(E, tabE);
98 1 : size_t j = random.randBin(tabCDF[i]);
99 1 : double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1);
100 1 : double eps = s_kin / 4. / E; // random background photon energy
101 :
102 : // Use approximation from A. Mastichiadis et al., Astroph. Journ. 300:178-189 (1986), eq. 30.
103 : // This approx is valid only for alpha >=100 where alpha = p0*eps*costheta - E0*eps
104 : // For our purposes, me << E0 --> p0~E0 --> alpha = E0*eps*(costheta - 1) >= 100
105 1 : double Epp = 5.7e-1 * pow(eps / mec2, -0.56) * pow(E / mec2, 0.44) * mec2;
106 :
107 1 : double f = Epp / E;
108 :
109 1 : if (haveElectrons) {
110 1 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
111 1 : if (random.rand() < pow(1 - f, thinning)) {
112 1 : double w = 1. / pow(1 - f, thinning);
113 1 : candidate->addSecondary(11, Epp / (1 + z), pos, w, interactionTag);
114 : }
115 1 : if (random.rand() < pow(f, thinning)) {
116 1 : double w = 1. / pow(f, thinning);
117 1 : candidate->addSecondary(-11, Epp / (1 + z), pos, w, interactionTag);
118 : }
119 : }
120 :
121 : // Update the primary particle energy.
122 : // This is done after adding the secondaries to correctly set the secondaries parent
123 1 : candidate->current.setEnergy((E - 2 * Epp) / (1. + z));
124 1 : }
125 :
126 1 : void EMTripletPairProduction::process(Candidate *candidate) const {
127 : // check if electron / positron
128 1 : int id = candidate->current.getId();
129 1 : if (abs(id) != 11)
130 1 : return;
131 :
132 : // scale the particle energy instead of background photons
133 1 : double z = candidate->getRedshift();
134 1 : double E = (1 + z) * candidate->current.getEnergy();
135 1 : Vector3d position = candidate->current.getPosition();
136 :
137 : // interaction rate
138 1 : double rate = getRate(E, position, z);
139 1 : if (rate < 0)
140 : return;
141 :
142 :
143 : // run this loop at least once to limit the step size
144 1 : double step = candidate->getCurrentStep();
145 1 : Random &random = Random::instance();
146 : do {
147 1 : double randDistance = -log(random.rand()) / rate;
148 : // check for interaction; if it doesn't occur, limit next step
149 1 : if (step < randDistance) {
150 1 : candidate->limitNextStep(limit / rate);
151 1 : return;
152 : }
153 0 : performInteraction(candidate);
154 0 : step -= randDistance;
155 0 : } while (step > 0.);
156 : }
157 :
158 1 : void EMTripletPairProduction::setInteractionTag(std::string tag) {
159 1 : interactionTag = tag;
160 1 : }
161 :
162 2 : std::string EMTripletPairProduction::getInteractionTag() const {
163 2 : return interactionTag;
164 : }
165 :
166 : } // namespace crpropa
|
