| 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 implementation file 29 // GEANT4 Class implementation file
30 // 30 //
31 // File name: G4GoudsmitSaundersonTable 31 // File name: G4GoudsmitSaundersonTable
32 // 32 //
33 // Author: Mihaly Novak / (Omrane Kadri 33 // Author: Mihaly Novak / (Omrane Kadri)
34 // 34 //
35 // Creation date: 20.02.2009 35 // Creation date: 20.02.2009
36 // 36 //
37 // Class description: 37 // Class description:
38 // Class to handle multiple scattering angul 38 // Class to handle multiple scattering angular distributions precomputed by
39 // using Kawrakow-Bielajew Goudsmit-Saunders 39 // using Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened
40 // Rutherford DCS for elastic scattering of 40 // Rutherford DCS for elastic scattering of electrons/positrons [1,2]. This
41 // class is used by G4GoudsmitSaundersonMscM 41 // class is used by G4GoudsmitSaundersonMscModel to sample the angular
42 // deflection of electrons/positrons after t 42 // deflection of electrons/positrons after travelling a given path.
43 // 43 //
44 // Modifications: 44 // Modifications:
45 // 04.03.2009 V.Ivanchenko cleanup and format 45 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
46 // 26.08.2009 O.Kadri: avoiding unuseful calcu 46 // 26.08.2009 O.Kadri: avoiding unuseful calculations and optimizing the root
47 // finding parameter error 47 // finding parameter error's within SampleTheta method
48 // 08.02.2010 O.Kadri: reduce delared variable 48 // 08.02.2010 O.Kadri: reduce delared variables; reduce error of finding root
49 // in secant method 49 // in secant method
50 // 26.03.2010 O.Kadri: minimum of used arrays 50 // 26.03.2010 O.Kadri: minimum of used arrays in computation within the dichotomie
51 // finding method the erro 51 // finding method the error was the lowest value of uvalues
52 // 12.05.2010 O.Kadri: changing of sqrt((b-a)* 52 // 12.05.2010 O.Kadri: changing of sqrt((b-a)*(b-a)) with fabs(b-a)
53 // 18.05.2015 M. Novak This class has been com 53 // 18.05.2015 M. Novak This class has been completely replaced (only the original
54 // class name was kept; class descr 54 // class name was kept; class description was also inserted):
55 // A new version of Kawrakow-Bielaj 55 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model
56 // based on the screened Rutherford 56 // based on the screened Rutherford DCS for elastic scattering of
57 // electrons/positrons has been int 57 // electrons/positrons has been introduced[1,2]. The corresponding MSC
58 // angular distributions over a 2D 58 // angular distributions over a 2D parameter grid have been recomputed
59 // and the CDFs are now stored in a 59 // and the CDFs are now stored in a variable transformed (smooth) form
60 // together with the corresponding 60 // together with the corresponding rational interpolation parameters.
61 // The new version is several times 61 // The new version is several times faster, more robust and accurate
62 // compared to the earlier version 62 // compared to the earlier version (G4GoudsmitSaundersonMscModel class
63 // that use these data has been als 63 // that use these data has been also completely replaced)
64 // 28.04.2017 M. Novak: New representation of 64 // 28.04.2017 M. Novak: New representation of the angular distribution data with
65 // significantly reduced data size. 65 // significantly reduced data size.
66 // 23.08.2017 M. Novak: Added funtionality to 66 // 23.08.2017 M. Novak: Added funtionality to handle Mott-correction to the
67 // base GS angular distributions an 67 // base GS angular distributions and some other factors (screening
68 // parameter, first and second mome 68 // parameter, first and second moments) when Mott-correction is
69 // activated in the GS-MSC model. 69 // activated in the GS-MSC model.
70 // 70 //
71 // References: 71 // References:
72 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-20 72 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-208
73 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19 73 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336
74 // 74 //
75 // ------------------------------------------- 75 // -----------------------------------------------------------------------------
76 76
77 #include "G4GoudsmitSaundersonTable.hh" 77 #include "G4GoudsmitSaundersonTable.hh"
78 78
79 79
80 #include "G4PhysicalConstants.hh" 80 #include "G4PhysicalConstants.hh"
81 #include "Randomize.hh" 81 #include "Randomize.hh"
82 #include "G4Log.hh" 82 #include "G4Log.hh"
83 #include "G4Exp.hh" 83 #include "G4Exp.hh"
84 84
85 #include "G4GSMottCorrection.hh" 85 #include "G4GSMottCorrection.hh"
86 #include "G4MaterialTable.hh" 86 #include "G4MaterialTable.hh"
87 #include "G4Material.hh" 87 #include "G4Material.hh"
88 #include "G4MaterialCutsCouple.hh" 88 #include "G4MaterialCutsCouple.hh"
89 #include "G4ProductionCutsTable.hh" 89 #include "G4ProductionCutsTable.hh"
90 #include "G4EmParameters.hh" 90 #include "G4EmParameters.hh"
91 91
92 #include "G4String.hh" 92 #include "G4String.hh"
93 93
94 #include <fstream> 94 #include <fstream>
95 #include <cstdlib> 95 #include <cstdlib>
96 #include <cmath> 96 #include <cmath>
97 97
98 #include <iostream> 98 #include <iostream>
99 #include <iomanip> 99 #include <iomanip>
100 100
101 // perecomputed GS angular distributions, base 101 // perecomputed GS angular distributions, based on the Screened-Rutherford DCS
102 // are the same for e- and e+ so make sure we 102 // are the same for e- and e+ so make sure we load them only onece
103 G4bool G4GoudsmitSaundersonTable::gIsInitialis 103 G4bool G4GoudsmitSaundersonTable::gIsInitialised = false;
104 // 104 //
105 std::vector<G4GoudsmitSaundersonTable::GSMSCAn 105 std::vector<G4GoudsmitSaundersonTable::GSMSCAngularDtr*> G4GoudsmitSaundersonTable::gGSMSCAngularDistributions1;
106 std::vector<G4GoudsmitSaundersonTable::GSMSCAn 106 std::vector<G4GoudsmitSaundersonTable::GSMSCAngularDtr*> G4GoudsmitSaundersonTable::gGSMSCAngularDistributions2;
107 // 107 //
108 std::vector<double> G4GoudsmitSaundersonTable: 108 std::vector<double> G4GoudsmitSaundersonTable::gMoliereBc;
109 std::vector<double> G4GoudsmitSaundersonTable: 109 std::vector<double> G4GoudsmitSaundersonTable::gMoliereXc2;
110 110
111 111
112 G4GoudsmitSaundersonTable::G4GoudsmitSaunderso 112 G4GoudsmitSaundersonTable::G4GoudsmitSaundersonTable(G4bool iselectron) {
113 fIsElectron = iselectron; 113 fIsElectron = iselectron;
114 // set initial values: final values will be 114 // set initial values: final values will be set in the Initialize method
115 fLogLambda0 = 0.; // wi 115 fLogLambda0 = 0.; // will be set properly at init.
116 fLogDeltaLambda = 0.; // wi 116 fLogDeltaLambda = 0.; // will be set properly at init.
117 fInvLogDeltaLambda = 0.; // wi 117 fInvLogDeltaLambda = 0.; // will be set properly at init.
118 fInvDeltaQ1 = 0.; // wi 118 fInvDeltaQ1 = 0.; // will be set properly at init.
119 fDeltaQ2 = 0.; // wi 119 fDeltaQ2 = 0.; // will be set properly at init.
120 fInvDeltaQ2 = 0.; // wi 120 fInvDeltaQ2 = 0.; // will be set properly at init.
121 // 121 //
122 fLowEnergyLimit = 0.1*CLHEP::keV; // w 122 fLowEnergyLimit = 0.1*CLHEP::keV; // will be set properly at init.
123 fHighEnergyLimit = 100.0*CLHEP::MeV; // w 123 fHighEnergyLimit = 100.0*CLHEP::MeV; // will be set properly at init.
124 // 124 //
125 fIsMottCorrection = false; // w 125 fIsMottCorrection = false; // will be set properly at init.
126 fIsPWACorrection = false; // w 126 fIsPWACorrection = false; // will be set properly at init.
127 fMottCorrection = nullptr; 127 fMottCorrection = nullptr;
128 // 128 //
129 fNumSPCEbinPerDec = 3; 129 fNumSPCEbinPerDec = 3;
130 } 130 }
131 131
132 G4GoudsmitSaundersonTable::~G4GoudsmitSaunders 132 G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable() {
133 for (std::size_t i=0; i<gGSMSCAngularDistrib 133 for (std::size_t i=0; i<gGSMSCAngularDistributions1.size(); ++i) {
134 if (gGSMSCAngularDistributions1[i]) { 134 if (gGSMSCAngularDistributions1[i]) {
135 delete [] gGSMSCAngularDistributions1[i] 135 delete [] gGSMSCAngularDistributions1[i]->fUValues;
136 delete [] gGSMSCAngularDistributions1[i] 136 delete [] gGSMSCAngularDistributions1[i]->fParamA;
137 delete [] gGSMSCAngularDistributions1[i] 137 delete [] gGSMSCAngularDistributions1[i]->fParamB;
138 delete gGSMSCAngularDistributions1[i]; 138 delete gGSMSCAngularDistributions1[i];
139 } 139 }
140 } 140 }
141 gGSMSCAngularDistributions1.clear(); 141 gGSMSCAngularDistributions1.clear();
142 for (std::size_t i=0; i<gGSMSCAngularDistrib 142 for (std::size_t i=0; i<gGSMSCAngularDistributions2.size(); ++i) {
143 if (gGSMSCAngularDistributions2[i]) { 143 if (gGSMSCAngularDistributions2[i]) {
144 delete [] gGSMSCAngularDistributions2[i] 144 delete [] gGSMSCAngularDistributions2[i]->fUValues;
145 delete [] gGSMSCAngularDistributions2[i] 145 delete [] gGSMSCAngularDistributions2[i]->fParamA;
146 delete [] gGSMSCAngularDistributions2[i] 146 delete [] gGSMSCAngularDistributions2[i]->fParamB;
147 delete gGSMSCAngularDistributions2[i]; 147 delete gGSMSCAngularDistributions2[i];
148 } 148 }
149 } 149 }
150 gGSMSCAngularDistributions2.clear(); 150 gGSMSCAngularDistributions2.clear();
151 if (fMottCorrection) { 151 if (fMottCorrection) {
152 delete fMottCorrection; 152 delete fMottCorrection;
153 fMottCorrection = nullptr; 153 fMottCorrection = nullptr;
154 } 154 }
155 // clear scp correction data 155 // clear scp correction data
156 for (std::size_t imc=0; imc<fSCPCPerMatCuts. 156 for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
157 if (fSCPCPerMatCuts[imc]) { 157 if (fSCPCPerMatCuts[imc]) {
158 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 158 fSCPCPerMatCuts[imc]->fVSCPC.clear();
159 delete fSCPCPerMatCuts[imc]; 159 delete fSCPCPerMatCuts[imc];
160 } 160 }
161 } 161 }
162 fSCPCPerMatCuts.clear(); 162 fSCPCPerMatCuts.clear();
163 // 163 //
164 gIsInitialised = false; 164 gIsInitialised = false;
165 } 165 }
166 166
167 void G4GoudsmitSaundersonTable::Initialise(G4d 167 void G4GoudsmitSaundersonTable::Initialise(G4double lownergylimit, G4double highenergylimit) {
168 fLowEnergyLimit = lownergylimit; 168 fLowEnergyLimit = lownergylimit;
169 fHighEnergyLimit = highenergylimit; 169 fHighEnergyLimit = highenergylimit;
170 G4double lLambdaMin = G4Log(gLAMBMIN); 170 G4double lLambdaMin = G4Log(gLAMBMIN);
171 G4double lLambdaMax = G4Log(gLAMBMAX); 171 G4double lLambdaMax = G4Log(gLAMBMAX);
172 fLogLambda0 = lLambdaMin; 172 fLogLambda0 = lLambdaMin;
173 fLogDeltaLambda = (lLambdaMax-lLambdaMin 173 fLogDeltaLambda = (lLambdaMax-lLambdaMin)/(gLAMBNUM-1.);
174 fInvLogDeltaLambda = 1./fLogDeltaLambda; 174 fInvLogDeltaLambda = 1./fLogDeltaLambda;
175 fInvDeltaQ1 = 1./((gQMAX1-gQMIN1)/(g 175 fInvDeltaQ1 = 1./((gQMAX1-gQMIN1)/(gQNUM1-1.));
176 fDeltaQ2 = (gQMAX2-gQMIN2)/(gQNUM 176 fDeltaQ2 = (gQMAX2-gQMIN2)/(gQNUM2-1.);
177 fInvDeltaQ2 = 1./fDeltaQ2; 177 fInvDeltaQ2 = 1./fDeltaQ2;
178 // load precomputed angular distributions an 178 // load precomputed angular distributions and set up several values used during the sampling
179 // these are particle independet => they go 179 // these are particle independet => they go to static container: load them only onece
180 if (!gIsInitialised) { 180 if (!gIsInitialised) {
181 // load pre-computed GS angular distributi 181 // load pre-computed GS angular distributions (computed based on Screened-Rutherford DCS)
182 LoadMSCData(); 182 LoadMSCData();
183 gIsInitialised = true; 183 gIsInitialised = true;
184 } 184 }
185 InitMoliereMSCParams(); 185 InitMoliereMSCParams();
186 // Mott-correction: particle(e- or e+) depen 186 // Mott-correction: particle(e- or e+) dependet so init them
187 if (fIsMottCorrection) { 187 if (fIsMottCorrection) {
188 if (!fMottCorrection) { 188 if (!fMottCorrection) {
189 fMottCorrection = new G4GSMottCorrection 189 fMottCorrection = new G4GSMottCorrection(fIsElectron);
190 } 190 }
191 fMottCorrection->Initialise(); 191 fMottCorrection->Initialise();
192 } 192 }
193 // init scattering power correction data; us 193 // init scattering power correction data; used only together with Mott-correction
194 // (Moliere's parameters must be initialised 194 // (Moliere's parameters must be initialised before)
195 if (fMottCorrection) { 195 if (fMottCorrection) {
196 InitSCPCorrection(); 196 InitSCPCorrection();
197 } 197 }
198 } 198 }
199 199
200 200
201 // samplig multiple scattering angles cos(thet 201 // samplig multiple scattering angles cos(theta) and sin(thata)
202 // - including no-scattering, single, "few" s 202 // - including no-scattering, single, "few" scattering cases as well
203 // - Mott-correction will be included if it w 203 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
204 // lambdaval : s/lambda_el 204 // lambdaval : s/lambda_el
205 // qval : s/lambda_el G1 205 // qval : s/lambda_el G1
206 // scra : screening parameter 206 // scra : screening parameter
207 // cost : will be the smapled cos(theta) 207 // cost : will be the smapled cos(theta)
208 // sint : will be the smapled sin(theta) 208 // sint : will be the smapled sin(theta)
209 // lekin : logarithm of the current kineti 209 // lekin : logarithm of the current kinetic energy
210 // beta2 : the corresponding beta square 210 // beta2 : the corresponding beta square
211 // matindx : index of the current material 211 // matindx : index of the current material
212 // returns true if it was msc 212 // returns true if it was msc
213 G4bool G4GoudsmitSaundersonTable::Sampling(G4d 213 G4bool G4GoudsmitSaundersonTable::Sampling(G4double lambdaval, G4double qval, G4double scra, G4double &cost,
214 G4d 214 G4double &sint, G4double lekin, G4double beta2, G4int matindx,
215 GSM 215 GSMSCAngularDtr **gsDtr, G4int &mcekini, G4int &mcdelti,
216 G4d 216 G4double &transfPar, G4bool isfirst) {
217 G4double rand0 = G4UniformRand(); 217 G4double rand0 = G4UniformRand();
218 G4double expn = G4Exp(-lambdaval); 218 G4double expn = G4Exp(-lambdaval);
219 // 219 //
220 // no scattering case 220 // no scattering case
221 if (rand0<expn) { 221 if (rand0<expn) {
222 cost = 1.0; 222 cost = 1.0;
223 sint = 0.0; 223 sint = 0.0;
224 return false; 224 return false;
225 } 225 }
226 // 226 //
227 // single scattering case : sample from the 227 // single scattering case : sample from the single scattering PDF
228 // - Mott-correction will be included if it 228 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
229 if (rand0<(1.+lambdaval)*expn) { 229 if (rand0<(1.+lambdaval)*expn) {
230 // cost is sampled in SingleScattering() 230 // cost is sampled in SingleScattering()
231 cost = SingleScattering(lambdaval, scra, l 231 cost = SingleScattering(lambdaval, scra, lekin, beta2, matindx);
232 // add protections 232 // add protections
233 if (cost<-1.0) cost = -1.0; 233 if (cost<-1.0) cost = -1.0;
234 if (cost>1.0) cost = 1.0; 234 if (cost>1.0) cost = 1.0;
235 // compute sin(theta) from the sampled cos 235 // compute sin(theta) from the sampled cos(theta)
236 G4double dum0 = 1.-cost; 236 G4double dum0 = 1.-cost;
237 sint = std::sqrt(dum0*(2.0-dum0)); 237 sint = std::sqrt(dum0*(2.0-dum0));
238 return false; 238 return false;
239 } 239 }
240 // 240 //
241 // handle this case: 241 // handle this case:
242 // -lambdaval < 1 i.e. mean #elastic ev 242 // -lambdaval < 1 i.e. mean #elastic events along the step is < 1 but
243 // the currently sampled case is not 0 243 // the currently sampled case is not 0 or 1 scattering. [Our minimal
244 // lambdaval (that we have precomputed 244 // lambdaval (that we have precomputed, transformed angular distributions
245 // stored in a form of equally probabe 245 // stored in a form of equally probabe intervalls together with rational
246 // interp. parameters) is 1.] 246 // interp. parameters) is 1.]
247 // -probability of having n elastic eve 247 // -probability of having n elastic events follows Poisson stat. with
248 // lambdaval parameter. 248 // lambdaval parameter.
249 // -the max. probability (when lambdava 249 // -the max. probability (when lambdaval=1) of having more than one
250 // elastic events is 0.2642411 and the 250 // elastic events is 0.2642411 and the prob of having 2,3,..,n elastic
251 // events decays rapidly with n. So se 251 // events decays rapidly with n. So set a max n to 10.
252 // -sampling of this cases is done in a 252 // -sampling of this cases is done in a one-by-one single elastic event way
253 // where the current #elastic event is 253 // where the current #elastic event is sampled from the Poisson distr.
254 if (lambdaval<1.0) { 254 if (lambdaval<1.0) {
255 G4double prob, cumprob; 255 G4double prob, cumprob;
256 prob = cumprob = expn; 256 prob = cumprob = expn;
257 G4double curcost,cursint; 257 G4double curcost,cursint;
258 // init cos(theta) and sin(theta) to the z 258 // init cos(theta) and sin(theta) to the zero scattering values
259 cost = 1.0; 259 cost = 1.0;
260 sint = 0.0; 260 sint = 0.0;
261 for (G4int iel=1; iel<10; ++iel) { 261 for (G4int iel=1; iel<10; ++iel) {
262 // prob of having iel scattering from Po 262 // prob of having iel scattering from Poisson
263 prob *= lambdaval/(G4double)iel; 263 prob *= lambdaval/(G4double)iel;
264 cumprob += prob; 264 cumprob += prob;
265 // 265 //
266 //sample cos(theta) from the singe scatt 266 //sample cos(theta) from the singe scattering pdf:
267 // - Mott-correction will be included if 267 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
268 curcost = SingleScattering(lambdav 268 curcost = SingleScattering(lambdaval, scra, lekin, beta2, matindx);
269 G4double dum0 = 1.-curcost; 269 G4double dum0 = 1.-curcost;
270 cursint = dum0*(2.0-dum0); // sin^ 270 cursint = dum0*(2.0-dum0); // sin^2(theta)
271 // 271 //
272 // if we got current deflection that is 272 // if we got current deflection that is not too small
273 // then update cos(theta) sin(theta) 273 // then update cos(theta) sin(theta)
274 if (cursint>1.0e-20) { 274 if (cursint>1.0e-20) {
275 cursint = std::sqrt(cursint); 275 cursint = std::sqrt(cursint);
276 G4double curphi = CLHEP::twopi*G4Unifo 276 G4double curphi = CLHEP::twopi*G4UniformRand();
277 cost = cost*curcost-sint*cu 277 cost = cost*curcost-sint*cursint*std::cos(curphi);
278 sint = std::sqrt(std::max(0 278 sint = std::sqrt(std::max(0.0, (1.0-cost)*(1.0+cost)));
279 } 279 }
280 // 280 //
281 // check if we have done enough scatteri 281 // check if we have done enough scattering i.e. sampling from the Poisson
282 if (rand0<cumprob) { 282 if (rand0<cumprob) {
283 return false; 283 return false;
284 } 284 }
285 } 285 }
286 // if reached the max iter i.e. 10 286 // if reached the max iter i.e. 10
287 return false; 287 return false;
288 } 288 }
289 // 289 //
290 // multiple scattering case with lambdavalue 290 // multiple scattering case with lambdavalue >= 1:
291 // - use the precomputed and transformed G 291 // - use the precomputed and transformed Goudsmit-Saunderson angular
292 // distributions to sample cos(theta) 292 // distributions to sample cos(theta)
293 // - Mott-correction will be included if i 293 // - Mott-correction will be included if it was requested by the user (i.e. if fIsMottCorrection=true)
294 cost = SampleCosTheta(lambdaval, qval, scra, 294 cost = SampleCosTheta(lambdaval, qval, scra, lekin, beta2, matindx, gsDtr, mcekini, mcdelti, transfPar, isfirst);
295 // add protections 295 // add protections
296 if (cost<-1.0) cost = -1.0; 296 if (cost<-1.0) cost = -1.0;
297 if (cost> 1.0) cost = 1.0; 297 if (cost> 1.0) cost = 1.0;
298 // compute cos(theta) and sin(theta) from th 298 // compute cos(theta) and sin(theta) from the sampled 1-cos(theta)
299 G4double dum0 = 1.0-cost; 299 G4double dum0 = 1.0-cost;
300 sint = std::sqrt(dum0*(2.0-dum0)); 300 sint = std::sqrt(dum0*(2.0-dum0));
301 // return true if it was msc 301 // return true if it was msc
302 return true; 302 return true;
303 } 303 }
304 304
305 305
306 G4double G4GoudsmitSaundersonTable::SampleCosT 306 G4double G4GoudsmitSaundersonTable::SampleCosTheta(G4double lambdaval, G4double qval, G4double scra,
307 307 G4double lekin, G4double beta2, G4int matindx,
308 308 GSMSCAngularDtr **gsDtr, G4int &mcekini,G4int &mcdelti,
309 309 G4double &transfPar, G4bool isfirst) {
310 G4double cost = 1.; 310 G4double cost = 1.;
311 // determine the base GS angular distributio 311 // determine the base GS angular distribution if it is the first call (when sub-step sampling is used)
312 if (isfirst) { 312 if (isfirst) {
313 *gsDtr = GetGSAngularDtr(scra, lambdaval, 313 *gsDtr = GetGSAngularDtr(scra, lambdaval, qval, transfPar);
314 } 314 }
315 // sample cost from the GS angular distribut 315 // sample cost from the GS angular distribution (computed based on Screened-Rutherford DCS)
316 cost = SampleGSSRCosTheta(*gsDtr, transfPar) 316 cost = SampleGSSRCosTheta(*gsDtr, transfPar);
317 // Mott-correction if it was requested by th 317 // Mott-correction if it was requested by the user
318 if (fIsMottCorrection && *gsDtr) { // no 318 if (fIsMottCorrection && *gsDtr) { // no Mott-correction in case of izotropic theta
319 static const G4int nlooplim = 1000; 319 static const G4int nlooplim = 1000;
320 G4int nloop = 0 ; // rejection loop 320 G4int nloop = 0 ; // rejection loop counter
321 // G4int ekindx = -1; // evaluate only 321 // G4int ekindx = -1; // evaluate only in the first call
322 // G4int deltindx = -1 ; // evaluate onl 322 // G4int deltindx = -1 ; // evaluate only in the first call
323 G4double val = fMottCorrection->GetMo 323 G4double val = fMottCorrection->GetMottRejectionValue(lekin, beta2, qval, cost, matindx, mcekini, mcdelti);
324 while (G4UniformRand()>val && ++nloop<nloo 324 while (G4UniformRand()>val && ++nloop<nlooplim) {
325 // sampling cos(theta) 325 // sampling cos(theta)
326 cost = SampleGSSRCosTheta(*gsDtr, transf 326 cost = SampleGSSRCosTheta(*gsDtr, transfPar);
327 val = fMottCorrection->GetMottRejection 327 val = fMottCorrection->GetMottRejectionValue(lekin, beta2, qval, cost, matindx, mcekini, mcdelti);
328 }; 328 };
329 } 329 }
330 return cost; 330 return cost;
331 } 331 }
332 332
333 333
334 // returns with cost sampled from the GS angul 334 // returns with cost sampled from the GS angular distribution computed based on Screened-Rutherford DCS
335 G4double G4GoudsmitSaundersonTable::SampleGSSR 335 G4double G4GoudsmitSaundersonTable::SampleGSSRCosTheta(const GSMSCAngularDtr *gsDtr, G4double transfpar) {
336 // check if isotropic theta (i.e. cost is un 336 // check if isotropic theta (i.e. cost is uniform on [-1:1])
337 if (!gsDtr) { 337 if (!gsDtr) {
338 return 1.-2.0*G4UniformRand(); 338 return 1.-2.0*G4UniformRand();
339 } 339 }
340 // 340 //
341 // sampling form the selected distribution 341 // sampling form the selected distribution
342 G4double ndatm1 = gsDtr->fNumData-1.; 342 G4double ndatm1 = gsDtr->fNumData-1.;
343 G4double delta = 1.0/ndatm1; 343 G4double delta = 1.0/ndatm1;
344 // determine lower cumulative bin inidex 344 // determine lower cumulative bin inidex
345 G4double rndm = G4UniformRand(); 345 G4double rndm = G4UniformRand();
346 G4int indxl = rndm*ndatm1; 346 G4int indxl = rndm*ndatm1;
347 G4double aval = rndm-indxl*delta; 347 G4double aval = rndm-indxl*delta;
348 G4double dum0 = delta*aval; 348 G4double dum0 = delta*aval;
349 349
350 G4double dum1 = (1.0+gsDtr->fParamA[indxl] 350 G4double dum1 = (1.0+gsDtr->fParamA[indxl]+gsDtr->fParamB[indxl])*dum0;
351 G4double dum2 = delta*delta + gsDtr->fPara 351 G4double dum2 = delta*delta + gsDtr->fParamA[indxl]*dum0 + gsDtr->fParamB[indxl]*aval*aval;
352 G4double sample = gsDtr->fUValues[indxl] + 352 G4double sample = gsDtr->fUValues[indxl] + dum1/dum2 *(gsDtr->fUValues[indxl+1]-gsDtr->fUValues[indxl]);
353 // transform back u to cos(theta) : 353 // transform back u to cos(theta) :
354 // this is the sampled cos(theta) = (2.0*par 354 // this is the sampled cos(theta) = (2.0*para*sample)/(1.0-sample+para)
355 return 1.-(2.0*transfpar*sample)/(1.0-sample 355 return 1.-(2.0*transfpar*sample)/(1.0-sample+transfpar);
356 } 356 }
357 357
358 358
359 // determine the GS angular distribution we ne 359 // determine the GS angular distribution we need to sample from: will set other things as well ...
360 G4GoudsmitSaundersonTable::GSMSCAngularDtr* G4 360 G4GoudsmitSaundersonTable::GSMSCAngularDtr* G4GoudsmitSaundersonTable::GetGSAngularDtr(G4double scra,
361 361 G4double &lambdaval, G4double &qval, G4double &transfpar) {
362 GSMSCAngularDtr *dtr = nullptr; 362 GSMSCAngularDtr *dtr = nullptr;
363 G4bool first = false; 363 G4bool first = false;
364 // isotropic cost above gQMAX2 (i.e. dtr sta 364 // isotropic cost above gQMAX2 (i.e. dtr stays nullptr)
365 if (qval<gQMAX2) { 365 if (qval<gQMAX2) {
366 G4int lamIndx = -1; // lambda value in 366 G4int lamIndx = -1; // lambda value index
367 G4int qIndx = -1; // lambda value in 367 G4int qIndx = -1; // lambda value index
368 // init to second grid Q values 368 // init to second grid Q values
369 G4int numQVal = gQNUM2; 369 G4int numQVal = gQNUM2;
370 G4double minQVal = gQMIN2; 370 G4double minQVal = gQMIN2;
371 G4double invDelQ = fInvDeltaQ2; 371 G4double invDelQ = fInvDeltaQ2;
372 G4double pIndxH = 0.; // probability of 372 G4double pIndxH = 0.; // probability of taking higher index
373 // check if first or second grid needs to 373 // check if first or second grid needs to be used
374 if (qval<gQMIN2) { // first grid 374 if (qval<gQMIN2) { // first grid
375 first = true; 375 first = true;
376 // protect against qval<gQMIN1 376 // protect against qval<gQMIN1
377 if (qval<gQMIN1) { 377 if (qval<gQMIN1) {
378 qval = gQMIN1; 378 qval = gQMIN1;
379 qIndx = 0; 379 qIndx = 0;
380 //pIndxH = 0.; 380 //pIndxH = 0.;
381 } 381 }
382 // set to first grid Q values 382 // set to first grid Q values
383 numQVal = gQNUM1; 383 numQVal = gQNUM1;
384 minQVal = gQMIN1; 384 minQVal = gQMIN1;
385 invDelQ = fInvDeltaQ1; 385 invDelQ = fInvDeltaQ1;
386 } 386 }
387 // make sure that lambda = s/lambda_el is 387 // make sure that lambda = s/lambda_el is in [gLAMBMIN,gLAMBMAX)
388 // lambda<gLAMBMIN=1 is already handeled b 388 // lambda<gLAMBMIN=1 is already handeled before so lambda>= gLAMBMIN for sure
389 if (lambdaval>=gLAMBMAX) { 389 if (lambdaval>=gLAMBMAX) {
390 lambdaval = gLAMBMAX-1.e-8; 390 lambdaval = gLAMBMAX-1.e-8;
391 lamIndx = gLAMBNUM-1; 391 lamIndx = gLAMBNUM-1;
392 } 392 }
393 G4double lLambda = G4Log(lambdaval); 393 G4double lLambda = G4Log(lambdaval);
394 // 394 //
395 // determine lower lambda (=s/lambda_el) i 395 // determine lower lambda (=s/lambda_el) index: linear interp. on log(lambda) scale
396 if (lamIndx<0) { 396 if (lamIndx<0) {
397 pIndxH = (lLambda-fLogLambda0)*fInvLogD 397 pIndxH = (lLambda-fLogLambda0)*fInvLogDeltaLambda;
398 lamIndx = (G4int)(pIndxH); // low 398 lamIndx = (G4int)(pIndxH); // lower index of the lambda bin
399 pIndxH = pIndxH-lamIndx; // proba 399 pIndxH = pIndxH-lamIndx; // probability of taking the higher index distribution
400 if (G4UniformRand()<pIndxH) { 400 if (G4UniformRand()<pIndxH) {
401 ++lamIndx; 401 ++lamIndx;
402 } 402 }
403 } 403 }
404 // 404 //
405 // determine lower Q (=s/lambda_el G1) ind 405 // determine lower Q (=s/lambda_el G1) index: linear interp. on Q
406 if (qIndx<0) { 406 if (qIndx<0) {
407 pIndxH = (qval-minQVal)*invDelQ; 407 pIndxH = (qval-minQVal)*invDelQ;
408 qIndx = (G4int)(pIndxH); // low 408 qIndx = (G4int)(pIndxH); // lower index of the Q bin
409 pIndxH = pIndxH-qIndx; 409 pIndxH = pIndxH-qIndx;
410 if (G4UniformRand()<pIndxH) { 410 if (G4UniformRand()<pIndxH) {
411 ++qIndx; 411 ++qIndx;
412 } 412 }
413 } 413 }
414 // set indx 414 // set indx
415 G4int indx = lamIndx*numQVal+qIndx; 415 G4int indx = lamIndx*numQVal+qIndx;
416 if (first) { 416 if (first) {
417 dtr = gGSMSCAngularDistributions1[indx]; 417 dtr = gGSMSCAngularDistributions1[indx];
418 } else { 418 } else {
419 dtr = gGSMSCAngularDistributions2[indx]; 419 dtr = gGSMSCAngularDistributions2[indx];
420 } 420 }
421 // dtr might be nullptr that indicates iso 421 // dtr might be nullptr that indicates isotropic cot distribution because:
422 // - if the selected lamIndx, qIndx corres 422 // - if the selected lamIndx, qIndx correspond to L(=s/lambda_el) and Q(=s/lambda_el G1) such that G1(=Q/L) > 1
423 // G1 should always be < 1 and if G1 is 423 // G1 should always be < 1 and if G1 is ~1 -> the dtr is isotropic (this can only happen in case of the 2. grid)
424 // 424 //
425 // compute the transformation parameter 425 // compute the transformation parameter
426 if (lambdaval>10.0) { 426 if (lambdaval>10.0) {
427 transfpar = 0.5*(-2.77164+lLambda*( 2.94 427 transfpar = 0.5*(-2.77164+lLambda*( 2.94874-lLambda*(0.1535754-lLambda*0.00552888) ));
428 } else { 428 } else {
429 transfpar = 0.5*(1.347+lLambda*(0.209364 429 transfpar = 0.5*(1.347+lLambda*(0.209364-lLambda*(0.45525-lLambda*(0.50142-lLambda*0.081234))));
430 } 430 }
431 transfpar *= (lambdaval+4.0)*scra; 431 transfpar *= (lambdaval+4.0)*scra;
432 } 432 }
433 // return with the selected GS angular distr 433 // return with the selected GS angular distribution that we need to sample cost from (if nullptr => isotropic cost)
434 return dtr; 434 return dtr;
435 } 435 }
436 436
437 437
438 void G4GoudsmitSaundersonTable::LoadMSCData() 438 void G4GoudsmitSaundersonTable::LoadMSCData() {
439 gGSMSCAngularDistributions1.resize(gLAMBNUM* 439 gGSMSCAngularDistributions1.resize(gLAMBNUM*gQNUM1,nullptr);
440 const G4String str1 = G4EmParameters::Instan 440 const G4String str1 = G4EmParameters::Instance()->GetDirLEDATA() + "/msc_GS/GSGrid_1/gsDistr_";
441 for (G4int il=0; il<gLAMBNUM; ++il) { 441 for (G4int il=0; il<gLAMBNUM; ++il) {
442 G4String fname = str1 + std::to_string(il) 442 G4String fname = str1 + std::to_string(il);
443 std::ifstream infile(fname,std::ios::in); 443 std::ifstream infile(fname,std::ios::in);
444 if (!infile.is_open()) { 444 if (!infile.is_open()) {
445 G4String msgc = "Cannot open file: " + f 445 G4String msgc = "Cannot open file: " + fname;
446 G4Exception("G4GoudsmitSaundersonTable:: 446 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
447 FatalException, msgc.c_str()); 447 FatalException, msgc.c_str());
448 return; 448 return;
449 } 449 }
450 for (G4int iq=0; iq<gQNUM1; ++iq) { 450 for (G4int iq=0; iq<gQNUM1; ++iq) {
451 auto gsd = new GSMSCAngularDtr(); 451 auto gsd = new GSMSCAngularDtr();
452 infile >> gsd->fNumData; 452 infile >> gsd->fNumData;
453 gsd->fUValues = new G4double[gsd->fNumDa 453 gsd->fUValues = new G4double[gsd->fNumData]();
454 gsd->fParamA = new G4double[gsd->fNumDa 454 gsd->fParamA = new G4double[gsd->fNumData]();
455 gsd->fParamB = new G4double[gsd->fNumDa 455 gsd->fParamB = new G4double[gsd->fNumData]();
456 G4double ddummy; 456 G4double ddummy;
457 infile >> ddummy; infile >> ddummy; 457 infile >> ddummy; infile >> ddummy;
458 for (G4int i=0; i<gsd->fNumData; ++i) { 458 for (G4int i=0; i<gsd->fNumData; ++i) {
459 infile >> gsd->fUValues[i]; 459 infile >> gsd->fUValues[i];
460 infile >> gsd->fParamA[i]; 460 infile >> gsd->fParamA[i];
461 infile >> gsd->fParamB[i]; 461 infile >> gsd->fParamB[i];
462 } 462 }
463 gGSMSCAngularDistributions1[il*gQNUM1+iq 463 gGSMSCAngularDistributions1[il*gQNUM1+iq] = gsd;
464 } 464 }
465 infile.close(); 465 infile.close();
466 } 466 }
467 // 467 //
468 // second grid 468 // second grid
469 gGSMSCAngularDistributions2.resize(gLAMBNUM* 469 gGSMSCAngularDistributions2.resize(gLAMBNUM*gQNUM2,nullptr);
470 const G4String str2 = G4EmParameters::Instan 470 const G4String str2 = G4EmParameters::Instance()->GetDirLEDATA() + "/msc_GS/GSGrid_2/gsDistr_";
471 for (G4int il=0; il<gLAMBNUM; ++il) { 471 for (G4int il=0; il<gLAMBNUM; ++il) {
472 G4String fname = str2 + std::to_string(il) 472 G4String fname = str2 + std::to_string(il);
473 std::ifstream infile(fname,std::ios::in); 473 std::ifstream infile(fname,std::ios::in);
474 if (!infile.is_open()) { 474 if (!infile.is_open()) {
475 G4String msgc = "Cannot open file: " + f 475 G4String msgc = "Cannot open file: " + fname;
476 G4Exception("G4GoudsmitSaundersonTable:: 476 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
477 FatalException, msgc.c_str()); 477 FatalException, msgc.c_str());
478 return; 478 return;
479 } 479 }
480 for (G4int iq=0; iq<gQNUM2; ++iq) { 480 for (G4int iq=0; iq<gQNUM2; ++iq) {
481 G4int numData; 481 G4int numData;
482 infile >> numData; 482 infile >> numData;
483 if (numData>1) { 483 if (numData>1) {
484 auto gsd = new GSMSCAngularDtr(); 484 auto gsd = new GSMSCAngularDtr();
485 gsd->fNumData = numData; 485 gsd->fNumData = numData;
486 gsd->fUValues = new G4double[gsd->fNum 486 gsd->fUValues = new G4double[gsd->fNumData]();
487 gsd->fParamA = new G4double[gsd->fNum 487 gsd->fParamA = new G4double[gsd->fNumData]();
488 gsd->fParamB = new G4double[gsd->fNum 488 gsd->fParamB = new G4double[gsd->fNumData]();
489 double ddummy; 489 double ddummy;
490 infile >> ddummy; infile >> ddummy; 490 infile >> ddummy; infile >> ddummy;
491 for (G4int i=0; i<gsd->fNumData; ++i) 491 for (G4int i=0; i<gsd->fNumData; ++i) {
492 infile >> gsd->fUValues[i]; 492 infile >> gsd->fUValues[i];
493 infile >> gsd->fParamA[i]; 493 infile >> gsd->fParamA[i];
494 infile >> gsd->fParamB[i]; 494 infile >> gsd->fParamB[i];
495 } 495 }
496 gGSMSCAngularDistributions2[il*gQNUM2+ 496 gGSMSCAngularDistributions2[il*gQNUM2+iq] = gsd;
497 } else { 497 } else {
498 gGSMSCAngularDistributions2[il*gQNUM2+ 498 gGSMSCAngularDistributions2[il*gQNUM2+iq] = nullptr;
499 } 499 }
500 } 500 }
501 infile.close(); 501 infile.close();
502 } 502 }
503 } 503 }
504 504
505 // samples cost in single scattering based on 505 // samples cost in single scattering based on Screened-Rutherford DCS
506 // (with Mott-correction if it was requested) 506 // (with Mott-correction if it was requested)
507 G4double G4GoudsmitSaundersonTable::SingleScat 507 G4double G4GoudsmitSaundersonTable::SingleScattering(G4double /*lambdaval*/, G4double scra,
508 508 G4double lekin, G4double beta2,
509 509 G4int matindx) {
510 G4double rand1 = G4UniformRand(); 510 G4double rand1 = G4UniformRand();
511 // sample cost from the Screened-Rutherford 511 // sample cost from the Screened-Rutherford DCS
512 G4double cost = 1.-2.0*scra*rand1/(1.0-rand 512 G4double cost = 1.-2.0*scra*rand1/(1.0-rand1+scra);
513 // Mott-correction if it was requested by th 513 // Mott-correction if it was requested by the user
514 if (fIsMottCorrection) { 514 if (fIsMottCorrection) {
515 static const G4int nlooplim = 1000; // re 515 static const G4int nlooplim = 1000; // rejection loop limit
516 G4int nloop = 0 ; // loop counter 516 G4int nloop = 0 ; // loop counter
517 G4int ekindx = -1 ; // evaluate only 517 G4int ekindx = -1 ; // evaluate only in the first call
518 G4int deltindx = 0 ; // single scatter 518 G4int deltindx = 0 ; // single scattering case
519 G4double q1 = 0.; // not used when 519 G4double q1 = 0.; // not used when deltindx = 0;
520 // computing Mott rejection function value 520 // computing Mott rejection function value
521 G4double val = fMottCorrection->GetMo 521 G4double val = fMottCorrection->GetMottRejectionValue(lekin, beta2, q1, cost,
522 522 matindx, ekindx, deltindx);
523 while (G4UniformRand()>val && ++nloop<nloo 523 while (G4UniformRand()>val && ++nloop<nlooplim) {
524 // sampling cos(theta) from the Screened 524 // sampling cos(theta) from the Screened-Rutherford DCS
525 rand1 = G4UniformRand(); 525 rand1 = G4UniformRand();
526 cost = 1.-2.0*scra*rand1/(1.0-rand1+scr 526 cost = 1.-2.0*scra*rand1/(1.0-rand1+scra);
527 // computing Mott rejection function val 527 // computing Mott rejection function value
528 val = fMottCorrection->GetMottRejectio 528 val = fMottCorrection->GetMottRejectionValue(lekin, beta2, q1, cost, matindx,
529 529 ekindx, deltindx);
530 }; 530 };
531 } 531 }
532 return cost; 532 return cost;
533 } 533 }
534 534
535 535
536 void G4GoudsmitSaundersonTable::GetMottCorrec 536 void G4GoudsmitSaundersonTable::GetMottCorrectionFactors(G4double logekin, G4double beta2,
537 537 G4int matindx, G4double &mcToScr,
538 538 G4double &mcToQ1, G4double &mcToG2PerG1) {
539 if (fIsMottCorrection) { 539 if (fIsMottCorrection) {
540 fMottCorrection->GetMottCorrectionFactors( 540 fMottCorrection->GetMottCorrectionFactors(logekin, beta2, matindx, mcToScr, mcToQ1, mcToG2PerG1);
541 } 541 }
542 } 542 }
543 543
544 544
545 // compute material dependent Moliere MSC para 545 // compute material dependent Moliere MSC parameters at initialisation
546 void G4GoudsmitSaundersonTable::InitMoliereMSC 546 void G4GoudsmitSaundersonTable::InitMoliereMSCParams() {
547 const G4double const1 = 7821.6; // [ 547 const G4double const1 = 7821.6; // [cm2/g]
548 const G4double const2 = 0.1569; // [ 548 const G4double const2 = 0.1569; // [cm2 MeV2 / g]
549 const G4double finstrc2 = 5.325135453E-5; / 549 const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
550 550
551 G4MaterialTable* theMaterialTable = G4Mater 551 G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
552 // get number of materials in the table 552 // get number of materials in the table
553 std::size_t numMaterials = theMaterialTable 553 std::size_t numMaterials = theMaterialTable->size();
554 // make sure that we have long enough vecto 554 // make sure that we have long enough vectors
555 if(gMoliereBc.size()<numMaterials) { 555 if(gMoliereBc.size()<numMaterials) {
556 gMoliereBc.resize(numMaterials); 556 gMoliereBc.resize(numMaterials);
557 gMoliereXc2.resize(numMaterials); 557 gMoliereXc2.resize(numMaterials);
558 } 558 }
559 G4double xi = 1.0; 559 G4double xi = 1.0;
560 G4int maxZ = 200; 560 G4int maxZ = 200;
561 if (fIsMottCorrection || fIsPWACorrection) 561 if (fIsMottCorrection || fIsPWACorrection) {
562 // xi = 1.0; <= always set to 1 from n 562 // xi = 1.0; <= always set to 1 from now on
563 maxZ = G4GSMottCorrection::GetMaxZet(); 563 maxZ = G4GSMottCorrection::GetMaxZet();
564 } 564 }
565 // 565 //
566 for (std::size_t imat=0; imat<numMaterials; 566 for (std::size_t imat=0; imat<numMaterials; ++imat) {
567 const G4Material* theMaterial = 567 const G4Material* theMaterial = (*theMaterialTable)[imat];
568 const G4ElementVector* theElemVect = 568 const G4ElementVector* theElemVect = theMaterial->GetElementVector();
569 const G4int numelems = 569 const G4int numelems = (G4int)theMaterial->GetNumberOfElements();
570 // 570 //
571 const G4double* theNbAtomsPerVolVe 571 const G4double* theNbAtomsPerVolVect = theMaterial->GetVecNbOfAtomsPerVolume();
572 G4double theTotNbAtomsPerVo 572 G4double theTotNbAtomsPerVol = theMaterial->GetTotNbOfAtomsPerVolume();
573 // 573 //
574 G4double zs = 0.0; 574 G4double zs = 0.0;
575 G4double zx = 0.0; 575 G4double zx = 0.0;
576 G4double ze = 0.0; 576 G4double ze = 0.0;
577 G4double sa = 0.0; 577 G4double sa = 0.0;
578 // 578 //
579 for(G4int ielem = 0; ielem < numelems; ie 579 for(G4int ielem = 0; ielem < numelems; ielem++) {
580 G4double zet = (*theElemVect)[ielem]->G 580 G4double zet = (*theElemVect)[ielem]->GetZ();
581 if (zet>maxZ) { 581 if (zet>maxZ) {
582 zet = (G4double)maxZ; 582 zet = (G4double)maxZ;
583 } 583 }
584 G4double iwa = (*theElemVect)[ielem]-> 584 G4double iwa = (*theElemVect)[ielem]->GetN();
585 G4double ipz = theNbAtomsPerVolVect[ie 585 G4double ipz = theNbAtomsPerVolVect[ielem]/theTotNbAtomsPerVol;
586 G4double dum = ipz*zet*(zet+xi); 586 G4double dum = ipz*zet*(zet+xi);
587 zs += dum; 587 zs += dum;
588 ze += dum*(-2.0/3.0)*G4Log(ze 588 ze += dum*(-2.0/3.0)*G4Log(zet);
589 zx += dum*G4Log(1.0+3.34*fins 589 zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet);
590 sa += ipz*iwa; 590 sa += ipz*iwa;
591 } 591 }
592 G4double density = theMaterial->GetDensit 592 G4double density = theMaterial->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
593 // 593 //
594 gMoliereBc[theMaterial->GetIndex()] = co 594 gMoliereBc[theMaterial->GetIndex()] = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
595 gMoliereXc2[theMaterial->GetIndex()] = co 595 gMoliereXc2[theMaterial->GetIndex()] = const2*density*zs/sa; // [MeV2/cm]
596 // change to Geant4 internal units of 1/l 596 // change to Geant4 internal units of 1/length and energ2/length
597 gMoliereBc[theMaterial->GetIndex()] *= 1 597 gMoliereBc[theMaterial->GetIndex()] *= 1.0/CLHEP::cm;
598 gMoliereXc2[theMaterial->GetIndex()] *= C 598 gMoliereXc2[theMaterial->GetIndex()] *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
599 } 599 }
600 } 600 }
601 601
602 602
603 // this method is temporary, will be removed/r 603 // this method is temporary, will be removed/replaced with a more effictien solution after 10.3.ref09
604 G4double G4GoudsmitSaundersonTable::ComputeSca 604 G4double G4GoudsmitSaundersonTable::ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut, G4double ekin) {
605 G4int imc = matcut->GetIndex(); 605 G4int imc = matcut->GetIndex();
606 G4double corFactor = 1.0; 606 G4double corFactor = 1.0;
607 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin< 607 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<=fSCPCPerMatCuts[imc]->fPrCut) {
608 return corFactor; 608 return corFactor;
609 } 609 }
610 // get the scattering power correction facto 610 // get the scattering power correction factor
611 G4double lekin = G4Log(ekin); 611 G4double lekin = G4Log(ekin);
612 G4double remaining = (lekin-fSCPCPerMatCuts 612 G4double remaining = (lekin-fSCPCPerMatCuts[imc]->fLEmin)*fSCPCPerMatCuts[imc]->fILDel;
613 std::size_t lindx = (std::size_t)remaining 613 std::size_t lindx = (std::size_t)remaining;
614 remaining -= lindx; 614 remaining -= lindx;
615 std::size_t imax = fSCPCPerMatCuts[imc]-> 615 std::size_t imax = fSCPCPerMatCuts[imc]->fVSCPC.size()-1;
616 if (lindx>=imax) { 616 if (lindx>=imax) {
617 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i 617 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[imax];
618 } else { 618 } else {
619 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l 619 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[lindx] + remaining*(fSCPCPerMatCuts[imc]->fVSCPC[lindx+1]-fSCPCPerMatCuts[imc]->fVSCPC[lindx]);
620 } 620 }
621 return corFactor; 621 return corFactor;
622 } 622 }
623 623
624 624
625 void G4GoudsmitSaundersonTable::InitSCPCorrect 625 void G4GoudsmitSaundersonTable::InitSCPCorrection() {
626 // get the material-cuts table 626 // get the material-cuts table
627 G4ProductionCutsTable *thePCTable = G4Produc 627 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
628 std::size_t numMatCuts = thePCTab 628 std::size_t numMatCuts = thePCTable->GetTableSize();
629 // clear container if any 629 // clear container if any
630 for (std::size_t imc=0; imc<fSCPCPerMatCuts. 630 for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
631 if (fSCPCPerMatCuts[imc]) { 631 if (fSCPCPerMatCuts[imc]) {
632 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 632 fSCPCPerMatCuts[imc]->fVSCPC.clear();
633 delete fSCPCPerMatCuts[imc]; 633 delete fSCPCPerMatCuts[imc];
634 fSCPCPerMatCuts[imc] = nullptr; 634 fSCPCPerMatCuts[imc] = nullptr;
635 } 635 }
636 } 636 }
637 // 637 //
638 // set size of the container and create the 638 // set size of the container and create the corresponding data structures
639 fSCPCPerMatCuts.resize(numMatCuts,nullptr); 639 fSCPCPerMatCuts.resize(numMatCuts,nullptr);
640 // loop over the material-cuts and create sc 640 // loop over the material-cuts and create scattering power correction data structure for each
641 for (G4int imc=0; imc<(G4int)numMatCuts; ++i 641 for (G4int imc=0; imc<(G4int)numMatCuts; ++imc) {
642 const G4MaterialCutsCouple *matCut = theP 642 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
643 // get e- production cut in the current ma 643 // get e- production cut in the current material-cuts in energy
644 G4double limit; 644 G4double limit;
645 G4double ecut; 645 G4double ecut;
646 if (fIsElectron) { 646 if (fIsElectron) {
647 ecut = (*(thePCTable->GetEnergyCutsVect 647 ecut = (*(thePCTable->GetEnergyCutsVector(idxG4ElectronCut)))[matCut->GetIndex()];
648 limit = 2.*ecut; 648 limit = 2.*ecut;
649 } else { 649 } else {
650 ecut = (*(thePCTable->GetEnergyCutsVect 650 ecut = (*(thePCTable->GetEnergyCutsVector(idxG4PositronCut)))[matCut->GetIndex()];
651 limit = ecut; 651 limit = ecut;
652 } 652 }
653 G4double min = std::max(limit,fLowEnergyLi 653 G4double min = std::max(limit,fLowEnergyLimit);
654 G4double max = fHighEnergyLimit; 654 G4double max = fHighEnergyLimit;
655 if (min>=max) { 655 if (min>=max) {
656 fSCPCPerMatCuts[imc] = new SCPCorrection 656 fSCPCPerMatCuts[imc] = new SCPCorrection();
657 fSCPCPerMatCuts[imc]->fIsUse = false; 657 fSCPCPerMatCuts[imc]->fIsUse = false;
658 fSCPCPerMatCuts[imc]->fPrCut = min; 658 fSCPCPerMatCuts[imc]->fPrCut = min;
659 continue; 659 continue;
660 } 660 }
661 G4int numEbins = fNumSPCEbinPerDec*G 661 G4int numEbins = fNumSPCEbinPerDec*G4lrint(std::log10(max/min));
662 numEbins = std::max(numEbins,3 662 numEbins = std::max(numEbins,3);
663 G4double lmin = G4Log(min); 663 G4double lmin = G4Log(min);
664 G4double ldel = G4Log(max/min)/(num 664 G4double ldel = G4Log(max/min)/(numEbins-1.0);
665 fSCPCPerMatCuts[imc] = new SCPCorrection() 665 fSCPCPerMatCuts[imc] = new SCPCorrection();
666 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi 666 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbins,1.0);
667 fSCPCPerMatCuts[imc]->fIsUse = true; 667 fSCPCPerMatCuts[imc]->fIsUse = true;
668 fSCPCPerMatCuts[imc]->fPrCut = min; 668 fSCPCPerMatCuts[imc]->fPrCut = min;
669 fSCPCPerMatCuts[imc]->fLEmin = lmin; 669 fSCPCPerMatCuts[imc]->fLEmin = lmin;
670 fSCPCPerMatCuts[imc]->fILDel = 1./ldel; 670 fSCPCPerMatCuts[imc]->fILDel = 1./ldel;
671 for (G4int ie=0; ie<numEbins; ++ie) { 671 for (G4int ie=0; ie<numEbins; ++ie) {
672 G4double ekin = G4Exp(lmin+ie*ldel); 672 G4double ekin = G4Exp(lmin+ie*ldel);
673 G4double scpCorr = 1.0; 673 G4double scpCorr = 1.0;
674 // compute correction factor: I.Kawrakow 674 // compute correction factor: I.Kawrakow NIMB 114(1996)307-326 (Eqs(32-37))
675 if (ie>0) { 675 if (ie>0) {
676 G4double tau = ekin/CLHEP::electr 676 G4double tau = ekin/CLHEP::electron_mass_c2;
677 G4double tauCut = ecut/CLHEP::electr 677 G4double tauCut = ecut/CLHEP::electron_mass_c2;
678 // Moliere's screening parameter 678 // Moliere's screening parameter
679 G4int matindx = (G4int)matCut->Get 679 G4int matindx = (G4int)matCut->GetMaterial()->GetIndex();
680 G4double A = GetMoliereXc2(mati 680 G4double A = GetMoliereXc2(matindx)/(4.0*tau*(tau+2.)*GetMoliereBc(matindx));
681 G4double gr = (1.+2.*A)*G4Log(1. 681 G4double gr = (1.+2.*A)*G4Log(1.+1./A)-2.;
682 G4double dum0 = (tau+2.)/(tau+1.); 682 G4double dum0 = (tau+2.)/(tau+1.);
683 G4double dum1 = tau+1.; 683 G4double dum1 = tau+1.;
684 G4double gm = G4Log(0.5*tau/tauC 684 G4double gm = G4Log(0.5*tau/tauCut) + (1.+dum0*dum0)*G4Log(2.*(tau-tauCut+2.)/(tau+4.))
685 - 0.25*(tau+2.)*( 685 - 0.25*(tau+2.)*( tau+2.+2.*(2.*tau+1.)/(dum1*dum1))*
686 G4Log((tau+4.)*( 686 G4Log((tau+4.)*(tau-tauCut)/tau/(tau-tauCut+2.))
687 + 0.5*(tau-2*tauCu 687 + 0.5*(tau-2*tauCut)*(tau+2.)*(1./(tau-tauCut)-1./(dum1*dum1));
688 if (gm<gr) { 688 if (gm<gr) {
689 gm = gm/gr; 689 gm = gm/gr;
690 } else { 690 } else {
691 gm = 1.; 691 gm = 1.;
692 } 692 }
693 G4double z0 = matCut->GetMaterial()-> 693 G4double z0 = matCut->GetMaterial()->GetIonisation()->GetZeffective();
694 scpCorr = 1.-gm*z0/(z0*(z0+1.)); 694 scpCorr = 1.-gm*z0/(z0*(z0+1.));
695 } 695 }
696 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo 696 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCorr;
697 } 697 }
698 } 698 }
699 } 699 }
700 700