| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class header file 29 // GEANT4 Class header file
30 // 30 //
31 // 31 //
32 // File name: G4eDPWAElasticDCS 32 // File name: G4eDPWAElasticDCS
33 // 33 //
34 // Author: Mihaly Novak 34 // Author: Mihaly Novak
35 // 35 //
36 // Creation date: 02.07.2020 36 // Creation date: 02.07.2020
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 39 //
40 // Class Description: 40 // Class Description:
41 // 41 //
42 // Contains numerical Differential Cross Secti 42 // Contains numerical Differential Cross Sections (DCS) for e-/e+ Coulomb
43 // scattering computed by Dirac Partial Wave A 43 // scattering computed by Dirac Partial Wave Analysis (DPWA) [1]:
44 // - electrostatic interaction, with a local e 44 // - electrostatic interaction, with a local exchange correction in the case of
45 // electrons (using Dirac-Fock e- densities; 45 // electrons (using Dirac-Fock e- densities; finite nuclear size with Fermi
46 // charge distribution; exchange potential w 46 // charge distribution; exchange potential with Furness and McCarthy for e-)[2]
47 // - correlation-polarization (projectiles cau 47 // - correlation-polarization (projectiles cause the polarization of the charge
48 // cloud of the target atom and the induced 48 // cloud of the target atom and the induced dipole moment acts back on the
49 // projectile) was accounted by using Local- 49 // projectile) was accounted by using Local-Density Approximation (LDA) [2]
50 // - absorption: not included since it's an in 50 // - absorption: not included since it's an inelastic channel [2] (the cor-
51 // responding excitations needs to be modell 51 // responding excitations needs to be modelled by a separate, independent,
52 // inelastic model). 52 // inelastic model).
53 // Using the above mentioned DPWA computation 53 // Using the above mentioned DPWA computation with a free atom approximation
54 // might lead to questionable results below fe 54 // might lead to questionable results below few hundred [eV] where possible
55 // solid state or bounding effects might start 55 // solid state or bounding effects might start to affect the potential.
56 // Nevertheless, the lower energy was set to 1 56 // Nevertheless, the lower energy was set to 10 eV in order to provide(at least)
57 // some model even at low energies (with this 57 // some model even at low energies (with this caution). The highest projectile
58 // kinetic energy is 100 [MeV]. 58 // kinetic energy is 100 [MeV].
59 // 59 //
60 // The class provides interface methods for el 60 // The class provides interface methods for elastic, first-, second-transport
61 // cross section computations as well as for s 61 // cross section computations as well as for sampling cosine of polar angular
62 // deflections. These interface methods are al 62 // deflections. These interface methods are also available for resricted cross
63 // section computations and angular deflection 63 // section computations and angular deflection sampling.
64 // 64 //
65 // References: 65 // References:
66 // 66 //
67 // [1] Salvat, F., Jablonski, A. and Powell, C 67 // [1] Salvat, F., Jablonski, A. and Powell, C.J., 2005. ELSEPA—Dirac partial-
68 // wave calculation of elastic scattering 68 // wave calculation of elastic scattering of electrons and positrons by
69 // atoms, positive ions and molecules. Com 69 // atoms, positive ions and molecules. Computer physics communications,
70 // 165(2), pp.157-190. 70 // 165(2), pp.157-190.
71 // [2] Salvat, F., 2003. Optical-model potenti 71 // [2] Salvat, F., 2003. Optical-model potential for electron and positron
72 // elastic scattering by atoms. Physical R 72 // elastic scattering by atoms. Physical Review A, 68(1), p.012708.
73 // [3] Benedito, E., Fernández-Varea, J.M. an 73 // [3] Benedito, E., Fernández-Varea, J.M. and Salvat, F.,2001. Mixed simulation
74 // of the multiple elastic scattering of e 74 // of the multiple elastic scattering of electrons and positrons using
75 // partial-wave differential cross-section 75 // partial-wave differential cross-sections. Nuclear Instruments and Methods
76 // in Physics Research Section B: Beam Int 76 // in Physics Research Section B: Beam Interactions with Materials and Atoms,
77 // 174(1-2), pp.91-110. 77 // 174(1-2), pp.91-110.
78 // 78 //
79 // ------------------------------------------- 79 // -------------------------------------------------------------------
80 80
81 #ifndef G4eDPWAElasticDCS_h 81 #ifndef G4eDPWAElasticDCS_h
82 #define G4eDPWAElasticDCS_h 1 82 #define G4eDPWAElasticDCS_h 1
83 83
84 84
85 #include <vector> 85 #include <vector>
86 #include <fstream> 86 #include <fstream>
87 #include <iomanip> 87 #include <iomanip>
88 #include <sstream> 88 #include <sstream>
89 89
90 #include "globals.hh" 90 #include "globals.hh"
91 #include "G4String.hh" 91 #include "G4String.hh"
92 #include "G4Physics2DVector.hh" 92 #include "G4Physics2DVector.hh"
93 93
94 94
95 #include "G4MaterialTable.hh" 95 #include "G4MaterialTable.hh"
96 #include "G4Material.hh" 96 #include "G4Material.hh"
97 #include "G4Element.hh" 97 #include "G4Element.hh"
98 #include "G4MaterialCutsCouple.hh" 98 #include "G4MaterialCutsCouple.hh"
99 #include "G4ProductionCutsTable.hh" 99 #include "G4ProductionCutsTable.hh"
100 100
101 101
102 #include "G4Log.hh" 102 #include "G4Log.hh"
103 #include "G4Exp.hh" 103 #include "G4Exp.hh"
104 104
105 105
106 class G4eDPWAElasticDCS { 106 class G4eDPWAElasticDCS {
107 107
108 public: 108 public:
109 109
110 // CTR: 110 // CTR:
111 // - iselectron : data for e- (for e+ othe 111 // - iselectron : data for e- (for e+ otherwise)
112 // - isrestricted : sampling of angular defl 112 // - isrestricted : sampling of angular deflection on restricted interavl is
113 // required (i.e. in case o 113 // required (i.e. in case of mixed-simulation models)
114 G4eDPWAElasticDCS(G4bool iselectron=true, G4 114 G4eDPWAElasticDCS(G4bool iselectron=true, G4bool isrestricted=false);
115 115
116 // DTR 116 // DTR
117 ~G4eDPWAElasticDCS(); 117 ~G4eDPWAElasticDCS();
118 118
119 // initialise for a given 'iz' atomic number 119 // initialise for a given 'iz' atomic number:
120 // - nothing happens if it has already been 120 // - nothing happens if it has already been initialised for that Z.
121 void InitialiseForZ(std::size_t iz); 121 void InitialiseForZ(std::size_t iz);
122 122
123 // Computes the elastic, first and second cr 123 // Computes the elastic, first and second cross sections for the given kinetic
124 // energy and target atom. 124 // energy and target atom.
125 // Cross sections are zero ff ekin is below/ 125 // Cross sections are zero ff ekin is below/above the kinetic energy grid
126 void ComputeCSPerAtom(G4int iz, G4double ek 126 void ComputeCSPerAtom(G4int iz, G4double ekin, G4double& elcs, G4double& tr1cs,
127 G4double& tr2cs, G4dou 127 G4double& tr2cs, G4double mumin=0.0, G4double mumax=1.0);
128 128
129 129
130 // samples cos(theta) i.e. cosine of the pol 130 // samples cos(theta) i.e. cosine of the polar angle of scattering in elastic
131 // interaction (Coulomb scattering) of the p 131 // interaction (Coulomb scattering) of the projectile (e- or e+ depending on
132 // fIsElectron) with kinetic energy of exp(' 132 // fIsElectron) with kinetic energy of exp('lekin'), target atom with atomic
133 // muber of 'iz'. See the 'SampleCosineTheta 133 // muber of 'iz'. See the 'SampleCosineThetaRestricted' for obtain samples on
134 // a restricted inteval. 134 // a restricted inteval.
135 G4double SampleCosineTheta(std::size_t iz, G 135 G4double SampleCosineTheta(std::size_t iz, G4double lekin, G4double r1,
136 G4double r2, G4do 136 G4double r2, G4double r3);
137 137
138 // samples cos(theta) i.e. cosine of the pol 138 // samples cos(theta) i.e. cosine of the polar angle of scattering in elastic
139 // interaction (Coulomb scattering) of the p 139 // interaction (Coulomb scattering) of the projectile (e- or e+ depending on
140 // fIsElectron) with kinetic energy of exp(' 140 // fIsElectron) with kinetic energy of exp('lekin'), target atom with atomic
141 // muber of 'iz'. 141 // muber of 'iz'.
142 // The cosine theta will be in the [costMin, 142 // The cosine theta will be in the [costMin, costMax] interval where costMin
143 // corresponds to a maximum allowed polar sc 143 // corresponds to a maximum allowed polar scattering angle thetaMax while
144 // costMin corresponds to minimum allowed po 144 // costMin corresponds to minimum allowed polar scatterin angle thetaMin.
145 // See the 'SampleCosineTheta' for obtain sa 145 // See the 'SampleCosineTheta' for obtain samples on the entire [-1,1] range.
146 G4double SampleCosineThetaRestricted(std::si 146 G4double SampleCosineThetaRestricted(std::size_t iz, G4double lekin,
147 G4doubl 147 G4double r1, G4double r2,
148 G4doubl 148 G4double costMax, G4double costMin);
149 149
150 // interpolate scattering power correction f 150 // interpolate scattering power correction form table buit at init.
151 G4double ComputeScatteringPowerCorrection(co 151 G4double ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut,
152 G4 152 G4double ekin);
153 153
154 // build scattering power correction table a 154 // build scattering power correction table at init.
155 void InitSCPCorrection(G4double lowEnerg 155 void InitSCPCorrection(G4double lowEnergyLimit, G4double highEnergyLimit);
156 156
157 157
158 private: 158 private:
159 159
160 // data structure to store one sampling tabl 160 // data structure to store one sampling table: combined Alias + RatIn
161 // NOTE: when Alias is used, sampling on a r 161 // NOTE: when Alias is used, sampling on a resctricted interval is not possible
162 // However, Alias makes possible faste 162 // However, Alias makes possible faster sampling. Alias is used in case
163 // of single scattering model while it 163 // of single scattering model while it's not used in case of mixed-model
164 // when restricted interval sampling i 164 // when restricted interval sampling is needed. This is controlled by
165 // the fIsRestrictedSamplingRequired f 165 // the fIsRestrictedSamplingRequired flag (false by default).
166 struct OneSamplingTable { 166 struct OneSamplingTable {
167 OneSamplingTable () = default; 167 OneSamplingTable () = default;
168 void SetSize(std::size_t nx, G4bool useAli 168 void SetSize(std::size_t nx, G4bool useAlias) {
169 fN = nx; 169 fN = nx;
170 // Alias 170 // Alias
171 if (useAlias) { 171 if (useAlias) {
172 fW.resize(nx); 172 fW.resize(nx);
173 fI.resize(nx); 173 fI.resize(nx);
174 } 174 }
175 // Ratin 175 // Ratin
176 fCum.resize(nx); 176 fCum.resize(nx);
177 fA.resize(nx); 177 fA.resize(nx);
178 fB.resize(nx); 178 fB.resize(nx);
179 } 179 }
180 180
181 // members 181 // members
182 std::size_t fN; // # 182 std::size_t fN; // # data points
183 G4double fScreenParA; // th 183 G4double fScreenParA; // the screening parameter
184 std::vector<G4double> fW; 184 std::vector<G4double> fW;
185 std::vector<G4double> fCum; 185 std::vector<G4double> fCum;
186 std::vector<G4double> fA; 186 std::vector<G4double> fA;
187 std::vector<G4double> fB; 187 std::vector<G4double> fB;
188 std::vector<G4int> fI; 188 std::vector<G4int> fI;
189 }; 189 };
190 190
191 191
192 // loads the kinetic energy and theta grids 192 // loads the kinetic energy and theta grids for the DCS data (first init step)
193 // should be called only by the master 193 // should be called only by the master
194 void LoadGrid(); 194 void LoadGrid();
195 195
196 // load DCS data for a given Z 196 // load DCS data for a given Z
197 void LoadDCSForZ(G4int iz); 197 void LoadDCSForZ(G4int iz);
198 198
199 // loads sampling table for the given Z over 199 // loads sampling table for the given Z over the enrgy grid
200 void BuildSmplingTableForZ(G4int iz); 200 void BuildSmplingTableForZ(G4int iz);
201 201
202 G4double SampleMu(std::size_t izet, std::siz 202 G4double SampleMu(std::size_t izet, std::size_t ie, G4double r1, G4double r2);
203 203
204 G4double FindCumValue(G4double u, const OneS 204 G4double FindCumValue(G4double u, const OneSamplingTable& stable,
205 const std::vector<G4do 205 const std::vector<G4double>& uvect);
206 206
207 // muMin and muMax : no checks on these 207 // muMin and muMax : no checks on these
208 G4double SampleMu(std::size_t izet, std::siz 208 G4double SampleMu(std::size_t izet, std::size_t ie, G4double r1, G4double muMin,
209 G4double muMax); 209 G4double muMax);
210 210
211 // set the DCS data directory path 211 // set the DCS data directory path
212 const G4String& FindDirectoryPath(); 212 const G4String& FindDirectoryPath();
213 213
214 // uncompress one data file into the input s 214 // uncompress one data file into the input string stream
215 void ReadCompressedFile(G4String fname, std: 215 void ReadCompressedFile(G4String fname, std::istringstream &iss);
216 216
217 // compute Molier material dependent paramet 217 // compute Molier material dependent parameters
218 void ComputeMParams(const G4Material* mat, G 218 void ComputeMParams(const G4Material* mat, G4double& theBc, G4double& theXc2);
219 219
220 220
221 // members 221 // members
222 private: 222 private:
223 223
224 // indicates if the object is for mixed-simu 224 // indicates if the object is for mixed-simulation (single scatterin otherwise)
225 G4bool fIsRestrict 225 G4bool fIsRestrictedSamplingRequired;
226 // indicates if the object is for e- (for e+ 226 // indicates if the object is for e- (for e+ otherwise)
227 G4bool fIsElectron 227 G4bool fIsElectron;
228 // indicates if the ekin, mu grids has alrea 228 // indicates if the ekin, mu grids has already been loaded (only once)
229 static G4bool gIsGridLoad 229 static G4bool gIsGridLoaded;
230 // data directory 230 // data directory
231 static G4String gDataDirect 231 static G4String gDataDirectory;
232 // max atomic number (Z) for which DCS has b 232 // max atomic number (Z) for which DCS has been computed (103)
233 static constexpr std::size_t gMaxZ = 103 233 static constexpr std::size_t gMaxZ = 103;
234 // energy and theta grid(s) relaed variables 234 // energy and theta grid(s) relaed variables: loaded from gridinfo by LoadGrid
235 static std::size_t gNumEnergie 235 static std::size_t gNumEnergies;
236 static std::size_t gIndxEnergy 236 static std::size_t gIndxEnergyLim;// the energy index just above 2 [keV]
237 static std::size_t gNumThetas1 237 static std::size_t gNumThetas1; // used for e- below 2 [keV]
238 static std::size_t gNumThetas2 238 static std::size_t gNumThetas2; // used for e+ and for e- bove 2 [keV]
239 static std::vector<G4double> gTheEnergie 239 static std::vector<G4double> gTheEnergies; // log-kinetic energy grid
240 static std::vector<G4double> gTheMus1; 240 static std::vector<G4double> gTheMus1; // mu(theta) = 0.5[1-cos(theta)]
241 static std::vector<G4double> gTheMus2; 241 static std::vector<G4double> gTheMus2;
242 static std::vector<G4double> gTheU1; 242 static std::vector<G4double> gTheU1; // u(mu; A'=0.01) = (A'+1)mu/(mu+A')
243 static std::vector<G4double> gTheU2; 243 static std::vector<G4double> gTheU2;
244 static G4double gLogMinEkin 244 static G4double gLogMinEkin; // log(gTheEnergies[0])
245 static G4double gInvDelLogE 245 static G4double gInvDelLogEkin;// 1./log(gTheEnergies[i+1]/gTheEnergies[i])
246 // abscissas and weights of an 8 point Gauss 246 // abscissas and weights of an 8 point Gauss-Legendre quadrature
247 // for numerical integration on [0,1] 247 // for numerical integration on [0,1]
248 static const G4double gXGL[8]; 248 static const G4double gXGL[8];
249 static const G4double gWGL[8]; 249 static const G4double gWGL[8];
250 // 250 //
251 std::vector<G4Physics2DVector*> fDCS; // 251 std::vector<G4Physics2DVector*> fDCS; // log(DCS) data per Z
252 std::vector<G4Physics2DVector*> fDCSLow; // 252 std::vector<G4Physics2DVector*> fDCSLow; // only for e- E < 2keV
253 // sampling tables: only one of the followin 253 // sampling tables: only one of the followings will be utilized
254 std::vector< std::vector<OneSamplingTable>* 254 std::vector< std::vector<OneSamplingTable>* > fSamplingTables;
255 // 255 //
256 // scattering power correction: to account s 256 // scattering power correction: to account sub-threshold inelastic deflections
257 const G4int fNumSPCEbin 257 const G4int fNumSPCEbinPerDec = 3;
258 struct SCPCorrection { 258 struct SCPCorrection {
259 G4bool fIsUse; // 259 G4bool fIsUse; //
260 G4double fPrCut; // sec. e- 260 G4double fPrCut; // sec. e- production cut energy
261 G4double fLEmin; // log min 261 G4double fLEmin; // log min energy
262 G4double fILDel; // inverse 262 G4double fILDel; // inverse log delta kinetic energy
263 std::vector<G4double> fVSCPC; // scatteri 263 std::vector<G4double> fVSCPC; // scattering power correction vector
264 }; 264 };
265 std::vector<SCPCorrection*> fSCPCPerMat 265 std::vector<SCPCorrection*> fSCPCPerMatCuts;
266 266
267 267
268 }; 268 };
269 269
270 #endif 270 #endif
271 271