Line data Source code
1 : #ifndef CRPROPA_PHOTOPIONPRODUCTION_H
2 : #define CRPROPA_PHOTOPIONPRODUCTION_H
3 :
4 : #include "crpropa/Module.h"
5 : #include "crpropa/PhotonBackground.h"
6 :
7 : #include <vector>
8 :
9 : namespace crpropa {
10 : /**
11 : * \addtogroup EnergyLosses
12 : * @{
13 : */
14 :
15 0 : struct SophiaEventOutput {
16 : int nParticles;
17 : std::vector<double> energy;
18 : std::vector<int> id;
19 : };
20 :
21 : /**
22 : @class PhotoPionProduction
23 : @brief Photo-pion interactions of nuclei with background photons.
24 : */
25 : class PhotoPionProduction: public Module {
26 :
27 : protected:
28 : ref_ptr<PhotonField> photonField;
29 : std::vector<double> tabLorentz; ///< Lorentz factor of nucleus
30 : std::vector<double> tabRedshifts; ///< redshifts (optional for haveRedshiftDependence)
31 : std::vector<double> tabProtonRate; ///< interaction rate in [1/m] for protons
32 : std::vector<double> tabNeutronRate; ///< interaction rate in [1/m] for neutrons
33 : double limit; ///< fraction of mean free path to limit the next step
34 : bool havePhotons;
35 : bool haveNeutrinos;
36 : bool haveElectrons;
37 : bool haveAntiNucleons;
38 : bool haveRedshiftDependence;
39 : std::string interactionTag = "PPP";
40 :
41 : // called by: sampleEps
42 : // - input: s [GeV^2]
43 : // - output: (s-p^2) * sigma_(nucleon/gamma) [GeV^2 * mubarn]
44 : double functs(double s, bool onProton) const;
45 :
46 : // called by: sampleEps, gaussInt
47 : // - input: photon energy eps [eV], Ein [GeV]
48 : // - output: probability to encounter photon of energy eps
49 : double probEps(double eps, bool onProton, double Ein, double z) const;
50 :
51 : /** called by: sampleEps
52 : @param onProton particle type: proton or neutron
53 : @param Ein energy of incoming nucleon
54 : - output: labframe energy [eV] of least energetic photon where PPP can occur
55 : */
56 : double epsMinInteraction(bool onProton, double Ein) const;
57 :
58 : /** called by: probEps, epsMinInteraction
59 : @param onProton particle type: proton or neutron
60 : @param Ein energy of incoming nucleon
61 : - output: hadron momentum [GeV/c]
62 : */
63 : double momentum(bool onProton, double Ein) const;
64 :
65 : // called by: functs
66 : // - input: photon energy [eV]
67 : // - output: crossection of nucleon-photon-interaction [mubarn]
68 : double crossection(double eps, bool onProton) const;
69 :
70 : // called by: crossection
71 : // - input: photon energy [eV], threshold [eV], max [eV], unknown [no unit]
72 : // - output: unknown [no unit]
73 : double Pl(double eps, double xth, double xMax, double alpha) const;
74 :
75 : // called by: crossection
76 : // - input: photon energy [eV], threshold [eV], unknown [eV]
77 : // - output: unknown [no unit]
78 : double Ef(double eps, double epsTh, double w) const;
79 :
80 : // called by: crossection
81 : // - input: cross section [µbarn], width [GeV], mass [GeV/c^2], rest frame photon energy [GeV]
82 : // - output: Breit-Wigner crossection of a resonance of width Gamma
83 : double breitwigner(double sigma0, double gamma, double DMM, double epsPrime, bool onProton) const;
84 :
85 : // called by: probEps, crossection, breitwigner, functs
86 : // - input: is proton [bool]
87 : // - output: mass [Gev/c^2]
88 : double mass(bool onProton) const;
89 :
90 : // - output: [GeV^2] head-on collision
91 : double sMin() const;
92 :
93 : bool sampleLog = true;
94 : double correctionFactor = 1.6; // increeses the maximum of the propability function
95 :
96 :
97 : public:
98 : /**
99 : * @brief pion production on a given target photon field
100 : *
101 : * @param photonField target photon field
102 : * @param photons if true, secondary photons are added to the simulation
103 : * @param neutrinos if true, secondary neutrinos are added to the simulation
104 : * @param electrons if true, secondary electrons are added to the simulation
105 : * @param antiNucleons if true, secondary anti nucleons are added to the simulation
106 : * @param limit fraction of the mean free path, to which the propagation step will be limited
107 : * @param haveRedshiftDependence use redshift dependent tabulated loss rates; if false, the redshift scaling of the photon field will be used
108 : */
109 : PhotoPionProduction(
110 : ref_ptr<PhotonField> photonField,
111 : bool photons = false,
112 : bool neutrinos = false,
113 : bool electrons = false,
114 : bool antiNucleons = false,
115 : double limit = 0.1,
116 : bool haveRedshiftDependence = false);
117 :
118 : // set the target photon field
119 : void setPhotonField(ref_ptr<PhotonField> photonField);
120 :
121 : // decide if secondary photons are added to the simulation
122 : void setHavePhotons(bool b);
123 :
124 : // decide if secondary neutrinos are added to the simulation
125 : void setHaveNeutrinos(bool b);
126 :
127 : // decide if secondary electrons are added to the simulation
128 : void setHaveElectrons(bool b);
129 :
130 : // decide if secondary anti nucleons are added to the simulation
131 : void setHaveAntiNucleons(bool b);
132 :
133 : // decide if redshift dependent tabulated loss rates are used
134 : void setHaveRedshiftDependence(bool b);
135 :
136 : /** Limit the propagation step to a fraction of the mean free path
137 : * @param limit fraction of the mean free path
138 : */
139 : void setLimit(double limit);
140 :
141 : /** set a custom interaction tag to trace back this interaction
142 : * @param tag string that will be added to the candidate and output
143 : */
144 : void setInteractionTag(std::string tag);
145 :
146 : void initRate(std::string filename);
147 :
148 : /** get the mean free path (MFP) for a single nucleon.
149 : * To get the MFP for the full nucleus the nucleonMFP has to be divided by by the nucleiModification factor
150 : * @param gamma Lorentz factor of the nucleon
151 : * @param z redshift
152 : * @param onProton true for protons, false for neutrons
153 : */
154 : double nucleonMFP(double gamma, double z, bool onProton) const;
155 :
156 : /** scaling factor for mean free path of the nucleus (converting the MFP of a single nucleon)
157 : *
158 : * @param A mass number of the nucleus
159 : * @param X charge number of the nucleus
160 : */
161 : double nucleiModification(int A, int X) const;
162 : void process(Candidate *candidate) const;
163 : void performInteraction(Candidate *candidate, bool onProton) const;
164 :
165 : /**
166 : Calculates the loss length E dx/dE in [m].
167 : This is not used in the simulation.
168 : @param id PDG particle id
169 : @param gamma Lorentz factor of particle
170 : @param z redshift
171 : */
172 : double lossLength(int id, double gamma, double z = 0);
173 :
174 : /**
175 : Direct SOPHIA interface.
176 : Output is an object SophiaEventOutput with two vectors "energy" and "id" each of length N (number of out-going particles).
177 : The i-th component of each vector corresponds to the same particle.
178 : This is not used in the simulation.
179 : @param onProton proton or neutron
180 : @param Ein energy of nucleon
181 : @param eps energy of target photon
182 : */
183 : SophiaEventOutput sophiaEvent(bool onProton, double Ein, double eps) const;
184 :
185 : /**
186 : SOPHIA's photon sampling method. Returns energy [J] of a photon of the photon field.
187 : @param onProton particle type: proton or neutron
188 : @param E energy of incoming nucleon [J]
189 : @param z redshift of incoming nucleon
190 : */
191 : double sampleEps(bool onProton, double E, double z) const;
192 :
193 : /** called by: sampleEps
194 : @param onProton particle type: proton or neutron
195 : @param Ein energy of incoming nucleon
196 : @param z redshift of incoming nucleon
197 : @param epsMin minimum photon energy of field
198 : @param epsMax maximum photon energy of field
199 : - output: maximum probability of all photons in field
200 : */
201 : double probEpsMax(bool onProton, double Ein, double z, double epsMin, double epsMax) const;
202 :
203 : // using log or lin spacing of photons in the range between epsMin and
204 : // epsMax for computing the maximum probability of photons in field
205 : void setSampleLog(bool log);
206 :
207 : // given the discrete steps to compute the maximum interaction probability pEpsMax
208 : // of photons in field, the real pEpsMax may lie between the descrete tested photon energies.
209 : // A correction factor can be set to increase pEpsMax by that factor
210 : void setCorrectionFactor(double factor);
211 :
212 : /** get functions for the parameters of the class PhotoPionProduction, similar to the set functions */
213 : ref_ptr<PhotonField> getPhotonField() const;
214 : bool getHavePhotons() const;
215 : bool getHaveNeutrinos() const;
216 : bool getHaveElectrons() const;
217 : bool getHaveAntiNucleons() const;
218 : bool getHaveRedshiftDependence() const;
219 : double getLimit() const;
220 : bool getSampleLog() const;
221 : double getCorrectionFactor() const;
222 : std::string getInteractionTag() const;
223 : };
224 : /** @}*/
225 :
226 : } // namespace crpropa
227 :
228 : #endif // CRPROPA_PHOTOPIONPRODUCTION_H
|
