| 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: G4GoudsmitSaundersonMscModel 31 // File name: G4GoudsmitSaundersonMscModel
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 // Modifications: 37 // Modifications:
38 // 04.03.2009 V.Ivanchenko cleanup and format 38 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as 39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as private doubles
40 // 18.05.2015 M. Novak provide PRELIMINARY ver 40 // 18.05.2015 M. Novak provide PRELIMINARY verison of the revised class
41 // This class has been revised and 41 // This class has been revised and updated, new methods added.
42 // A new version of Kawrakow-Bielaj 42 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model
43 // based on the screened Rutherford 43 // based on the screened Rutherford DCS for elastic scattering of
44 // electrons/positrons has been int 44 // electrons/positrons has been introduced[1,2]. The corresponding MSC
45 // angular distributions over a 2D 45 // angular distributions over a 2D parameter grid have been recomputed
46 // and the CDFs are now stored in a 46 // and the CDFs are now stored in a variable transformed (smooth) form[2,3]
47 // together with the corresponding 47 // together with the corresponding rational interpolation parameters.
48 // These angular distributions are 48 // These angular distributions are handled by the new
49 // G4GoudsmitSaundersonTable class 49 // G4GoudsmitSaundersonTable class that is responsible to sample if
50 // it was no, single, few or multip 50 // it was no, single, few or multiple scattering case and delivers the
51 // angular deflection (i.e. cos(the 51 // angular deflection (i.e. cos(theta) and sin(theta)).
52 // Two screening options are provid 52 // Two screening options are provided:
53 // - if fIsUsePWATotalXsecData=TRU 53 // - if fIsUsePWATotalXsecData=TRUE i.e. SetOptionPWAScreening(TRUE)
54 // was called before initialisat 54 // was called before initialisation: screening parameter value A is
55 // determined such that the firs 55 // determined such that the first transport coefficient G1(A)
56 // computed according to the scr 56 // computed according to the screened Rutherford DCS for elastic
57 // scattering will reproduce the 57 // scattering will reproduce the one computed from the PWA elastic
58 // and first transport mean free 58 // and first transport mean free paths[4].
59 // - if fIsUsePWATotalXsecData=FAL 59 // - if fIsUsePWATotalXsecData=FALSE i.e. default value or
60 // SetOptionPWAScreening(FALSE) 60 // SetOptionPWAScreening(FALSE) was called before initialisation:
61 // screening parameter value A i 61 // screening parameter value A is computed according to Moliere's
62 // formula (by using material de 62 // formula (by using material dependent parameters \chi_cc2 and b_c
63 // precomputed for each material 63 // precomputed for each material used at initialization in
64 // G4GoudsmitSaundersonTable) [3 64 // G4GoudsmitSaundersonTable) [3]
65 // Elastic and first trasport mean 65 // Elastic and first trasport mean free paths are used consistently.
66 // The new version is self-consiste 66 // The new version is self-consistent, several times faster, more
67 // robust and accurate compared to 67 // robust and accurate compared to the earlier version.
68 // Spin effects as well as a more a 68 // Spin effects as well as a more accurate energy loss correction and
69 // computations of Lewis moments wi 69 // computations of Lewis moments will be implemented later on.
70 // 02.09.2015 M. Novak: first version of new s 70 // 02.09.2015 M. Novak: first version of new step limit is provided.
71 // fUseSafetyPlus corresponds to Ur 71 // fUseSafetyPlus corresponds to Urban fUseSafety (default)
72 // fUseDistanceToBoundary correspon 72 // fUseDistanceToBoundary corresponds to Urban fUseDistanceToBoundary
73 // fUseSafety corresponds to EGSnr 73 // fUseSafety corresponds to EGSnrc error-free stepping algorithm
74 // Range factor can be significantl 74 // Range factor can be significantly higher at each case than in Urban.
75 // 23.08.2017 M. Novak: added corrections to a 75 // 23.08.2017 M. Novak: added corrections to account spin effects (Mott-correction).
76 // It can be activated by setting t 76 // It can be activated by setting the fIsMottCorrection flag to be true
77 // before initialization using the 77 // before initialization using the SetOptionMottCorrection() public method.
78 // The fMottCorrection member is re 78 // The fMottCorrection member is responsible to handle pre-computed Mott
79 // correction (rejection) functions 79 // correction (rejection) functions obtained by numerically computing
80 // Goudsmit-Saunderson agnular dist 80 // Goudsmit-Saunderson agnular distributions based on a DCS accounting spin
81 // effects and screening correction 81 // effects and screening corrections. The DCS used to compute the accurate
82 // GS angular distributions is: DCS 82 // GS angular distributions is: DCS_{cor} = DCS_{SR}x[ DCS_{R}/DCS_{Mott}] where :
83 // # DCS_{SR} is the relativisti 83 // # DCS_{SR} is the relativistic Screened-Rutherford DCS (first Born approximate
84 // solution of the Klein-Gordo 84 // solution of the Klein-Gordon i.e. relativistic Schrodinger equation =>
85 // scattering of spinless e- o 85 // scattering of spinless e- on exponentially screened Coulomb potential)
86 // note: the default (without 86 // note: the default (without using Mott-correction) GS angular distributions
87 // are based on this DCS_{SR} 87 // are based on this DCS_{SR} with Moliere's screening parameter!
88 // # DCS_{R} is the Rutherford D 88 // # DCS_{R} is the Rutherford DCS which is the same as above but without
89 // screening 89 // screening
90 // # DCS_{Mott} is the Mott DCS 90 // # DCS_{Mott} is the Mott DCS i.e. solution of the Dirac equation with a bare
91 // Coulomb potential i.e. scat 91 // Coulomb potential i.e. scattering of particles with spin (e- or e+) on a
92 // point-like unscreened Coulo 92 // point-like unscreened Coulomb potential
93 // # moreover, the screening par 93 // # moreover, the screening parameter of the DCS_{cor} was determined such that
94 // the DCS_{cor} with this cor 94 // the DCS_{cor} with this corrected screening parameter reproduce the first
95 // transport cross sections ob 95 // transport cross sections obtained from the corresponding most accurate DCS
96 // (i.e. from elsepa [4]) 96 // (i.e. from elsepa [4])
97 // Unlike the default GS, the Mott- 97 // Unlike the default GS, the Mott-corrected angular distributions are particle type
98 // (different for e- and e+ <= the 98 // (different for e- and e+ <= the DCS_{Mott} and the screening correction) and target
99 // (Z and material) dependent. 99 // (Z and material) dependent.
100 // 27.10.2017 M. Novak: 100 // 27.10.2017 M. Novak:
101 // - Mott-correction flag is set no 101 // - Mott-correction flag is set now through the G4EmParameters
102 // - new form of PWA correction to 102 // - new form of PWA correction to integrated quantities and screening (default)
103 // - changed step limit flag conven 103 // - changed step limit flag conventions:
104 // # fUseSafety corresponds to U 104 // # fUseSafety corresponds to Urban's fUseSafety
105 // # fUseDistanceToBoundary corr 105 // # fUseDistanceToBoundary corresponds to Urban's fUseDistanceToBoundary
106 // # fUseSafetyPlus corresponds 106 // # fUseSafetyPlus corresponds to the error-free stepping algorithm
107 // 02.02.2018 M. Novak: implemented CrossSecti 107 // 02.02.2018 M. Novak: implemented CrossSectionPerVolume interface method (used only for testing)
108 // 108 //
109 // Class description: 109 // Class description:
110 // Kawrakow-Bielajew Goudsmit-Saunderson MSC 110 // Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened Rutherford DCS
111 // for elastic scattering of e-/e+. Option, 111 // for elastic scattering of e-/e+. Option, to include (Mott) correction (see above), is
112 // also available now (SetOptionMottCorrecti 112 // also available now (SetOptionMottCorrection(true)). An EGSnrc like error-free stepping
113 // algorithm (UseSafety) is available beyond 113 // algorithm (UseSafety) is available beyond the usual Geant4 step limitation algorithms
114 // and true to geomerty and geometry to true 114 // and true to geomerty and geometry to true step length computations that were adopted
115 // from the Urban model[5]. The most accurat 115 // from the Urban model[5]. The most accurate setting: error-free stepping i.e. the
116 // UseSafetyPlus MSC step limit with Mott-co 116 // UseSafetyPlus MSC step limit with Mott-correction (SetOptionMottCorrection(true)). Both
117 // are expected to be set through the G4EmPa 117 // are expected to be set through the G4EmParameters singleton before initialisation:
118 // # G4EmParameters::Instance()->SetMscStep 118 // # G4EmParameters::Instance()->SetMscStepLimitType(fUseSafetyPlus);
119 // # G4EmParameters::Instance()->SetUseMott 119 // # G4EmParameters::Instance()->SetUseMottCorrection(true);
120 // 120 //
121 // 121 //
122 // References: 122 // References:
123 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208 123 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208
124 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19 124 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336
125 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Ro 125 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Rogers, F.Tessier,B.R.B.Walters, NRCC
126 // Report PIRS-701 (2013) 126 // Report PIRS-701 (2013)
127 // [4] F.Salvat, A.Jablonski, C.J. Powell, C 127 // [4] F.Salvat, A.Jablonski, C.J. Powell, CPC 165(2005) 157-190
128 // [5] L.Urban, Preprint CERN-OPEN-2006-077 128 // [5] L.Urban, Preprint CERN-OPEN-2006-077 (2006)
129 // 129 //
130 // ------------------------------------------- 130 // -----------------------------------------------------------------------------
131 131
132 132
133 #include "G4GoudsmitSaundersonMscModel.hh" 133 #include "G4GoudsmitSaundersonMscModel.hh"
134 134
135 #include "G4GoudsmitSaundersonTable.hh" 135 #include "G4GoudsmitSaundersonTable.hh"
136 #include "G4GSPWACorrections.hh" 136 #include "G4GSPWACorrections.hh"
137 137
138 #include "G4PhysicalConstants.hh" 138 #include "G4PhysicalConstants.hh"
139 #include "G4SystemOfUnits.hh" 139 #include "G4SystemOfUnits.hh"
140 140
141 #include "G4ParticleChangeForMSC.hh" 141 #include "G4ParticleChangeForMSC.hh"
142 #include "G4DynamicParticle.hh" 142 #include "G4DynamicParticle.hh"
143 #include "G4Electron.hh" 143 #include "G4Electron.hh"
144 #include "G4Positron.hh" 144 #include "G4Positron.hh"
145 145
146 #include "G4LossTableManager.hh" 146 #include "G4LossTableManager.hh"
147 #include "G4EmParameters.hh" 147 #include "G4EmParameters.hh"
148 #include "G4Track.hh" 148 #include "G4Track.hh"
149 #include "G4PhysicsTable.hh" 149 #include "G4PhysicsTable.hh"
150 #include "Randomize.hh" 150 #include "Randomize.hh"
151 #include "G4Log.hh" 151 #include "G4Log.hh"
152 #include "G4Exp.hh" 152 #include "G4Exp.hh"
153 #include "G4Pow.hh" 153 #include "G4Pow.hh"
154 #include <fstream> 154 #include <fstream>
155 155
156 156
157 // set accurate energy loss and dispalcement s 157 // set accurate energy loss and dispalcement sampling to be always on now
158 G4bool G4GoudsmitSaundersonMscModel::gIsUseAcc 158 G4bool G4GoudsmitSaundersonMscModel::gIsUseAccurate = true;
159 // set the usual optimization to be always act 159 // set the usual optimization to be always active now
160 G4bool G4GoudsmitSaundersonMscModel::gIsOptimi 160 G4bool G4GoudsmitSaundersonMscModel::gIsOptimizationOn = true;
161 161
162 162
163 G4GoudsmitSaundersonMscModel::G4GoudsmitSaunde 163 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
164 : G4VMscModel(nam) { 164 : G4VMscModel(nam) {
165 charge = 0; 165 charge = 0;
166 currentMaterialIndex = -1; 166 currentMaterialIndex = -1;
167 // 167 //
168 fr = 0.1; 168 fr = 0.1;
169 rangeinit = 1.e+21; 169 rangeinit = 1.e+21;
170 geombig = 1.e+50*mm; 170 geombig = 1.e+50*mm;
171 geomlimit = geombig; 171 geomlimit = geombig;
172 tgeom = geombig; 172 tgeom = geombig;
173 tlimit = 1.e+10*mm; 173 tlimit = 1.e+10*mm;
174 presafety = 0.*mm; 174 presafety = 0.*mm;
175 // 175 //
176 particle = nullptr; << 176 particle = 0;
177 theManager = G4LossTableManager: 177 theManager = G4LossTableManager::Instance();
178 firstStep = true; 178 firstStep = true;
179 currentKinEnergy = 0.0; 179 currentKinEnergy = 0.0;
180 currentRange = 0.0; 180 currentRange = 0.0;
181 // 181 //
182 tlimitminfix2 = 1.*nm; 182 tlimitminfix2 = 1.*nm;
183 tausmall = 1.e-16; 183 tausmall = 1.e-16;
184 mass = electron_mass_c2; 184 mass = electron_mass_c2;
185 taulim = 1.e-6; 185 taulim = 1.e-6;
186 // 186 //
187 currentCouple = nullptr; 187 currentCouple = nullptr;
188 fParticleChange = nullptr; 188 fParticleChange = nullptr;
189 // 189 //
190 fZeff = 1.; 190 fZeff = 1.;
191 // 191 //
192 par1 = 0.; 192 par1 = 0.;
193 par2 = 0.; 193 par2 = 0.;
194 par3 = 0.; 194 par3 = 0.;
195 // 195 //
196 // Moliere screeing parameter will be used a 196 // Moliere screeing parameter will be used and (by default) corrections are
197 // appalied to the integrated quantities (sc 197 // appalied to the integrated quantities (screeing parameter, elastic mfp, first
198 // and second moments) derived from the corr 198 // and second moments) derived from the corresponding PWA quantities
199 // this PWA correction is ignored if Mott-co 199 // this PWA correction is ignored if Mott-correction is set to true because
200 // Mott-correction contains all these correc 200 // Mott-correction contains all these corrections as well
201 fIsUsePWACorrection = true; 201 fIsUsePWACorrection = true;
202 // 202 //
203 fIsUseMottCorrection = false; 203 fIsUseMottCorrection = false;
204 // 204 //
205 fLambda0 = 0.0; // elastic mea 205 fLambda0 = 0.0; // elastic mean free path
206 fLambda1 = 0.0; // first trans 206 fLambda1 = 0.0; // first transport mean free path
207 fScrA = 0.0; // screening p 207 fScrA = 0.0; // screening parameter
208 fG1 = 0.0; // first trans 208 fG1 = 0.0; // first transport coef.
209 // 209 //
210 fMCtoScrA = 1.0; 210 fMCtoScrA = 1.0;
211 fMCtoQ1 = 1.0; 211 fMCtoQ1 = 1.0;
212 fMCtoG2PerG1 = 1.0; 212 fMCtoG2PerG1 = 1.0;
213 // 213 //
214 fTheTrueStepLenght = 0.; 214 fTheTrueStepLenght = 0.;
215 fTheTransportDistance = 0.; 215 fTheTransportDistance = 0.;
216 fTheZPathLenght = 0.; 216 fTheZPathLenght = 0.;
217 // 217 //
218 fTheDisplacementVector.set(0.,0.,0.); 218 fTheDisplacementVector.set(0.,0.,0.);
219 fTheNewDirection.set(0.,0.,1.); 219 fTheNewDirection.set(0.,0.,1.);
220 // 220 //
221 fIsEverythingWasDone = false; 221 fIsEverythingWasDone = false;
222 fIsMultipleSacettring = false; 222 fIsMultipleSacettring = false;
223 fIsSingleScattering = false; 223 fIsSingleScattering = false;
224 fIsEndedUpOnBoundary = false; 224 fIsEndedUpOnBoundary = false;
225 fIsNoScatteringInMSC = false; 225 fIsNoScatteringInMSC = false;
226 fIsNoDisplace = false; 226 fIsNoDisplace = false;
227 fIsInsideSkin = false; 227 fIsInsideSkin = false;
228 fIsWasOnBoundary = false; 228 fIsWasOnBoundary = false;
229 fIsFirstRealStep = false; 229 fIsFirstRealStep = false;
230 rndmEngineMod = G4Random::getTheEng 230 rndmEngineMod = G4Random::getTheEngine();
231 // 231 //
232 fGSTable = nullptr; 232 fGSTable = nullptr;
233 fPWACorrection = nullptr; 233 fPWACorrection = nullptr;
234 } 234 }
235 235
236 236
237 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaund 237 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel() {
238 if (IsMaster()) { 238 if (IsMaster()) {
239 if (fGSTable) { 239 if (fGSTable) {
240 delete fGSTable; 240 delete fGSTable;
241 fGSTable = nullptr; 241 fGSTable = nullptr;
242 } 242 }
243 if (fPWACorrection) { 243 if (fPWACorrection) {
244 delete fPWACorrection; 244 delete fPWACorrection;
245 fPWACorrection = nullptr; 245 fPWACorrection = nullptr;
246 } 246 }
247 } 247 }
248 } 248 }
249 249
250 250
251 void G4GoudsmitSaundersonMscModel::Initialise( 251 void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p, const G4DataVector&) {
252 SetParticle(p); 252 SetParticle(p);
253 InitialiseParameters(p); 253 InitialiseParameters(p);
254 // -create GoudsmitSaundersonTable and init 254 // -create GoudsmitSaundersonTable and init its Mott-correction member if
255 // Mott-correction was required 255 // Mott-correction was required
256 if (IsMaster()) { 256 if (IsMaster()) {
257 // get the Mott-correction flag from EmPar 257 // get the Mott-correction flag from EmParameters
258 if (G4EmParameters::Instance()->UseMottCor 258 if (G4EmParameters::Instance()->UseMottCorrection()) {
259 fIsUseMottCorrection = true; 259 fIsUseMottCorrection = true;
260 } 260 }
261 // Mott-correction includes other way of P 261 // Mott-correction includes other way of PWA x-section corrections so deactivate it even if it was true
262 // when Mott-correction is activated by th 262 // when Mott-correction is activated by the user
263 if (fIsUseMottCorrection) { 263 if (fIsUseMottCorrection) {
264 fIsUsePWACorrection = false; 264 fIsUsePWACorrection = false;
265 } 265 }
266 // clear GS-table 266 // clear GS-table
267 if (fGSTable) { 267 if (fGSTable) {
268 delete fGSTable; 268 delete fGSTable;
269 fGSTable = nullptr; 269 fGSTable = nullptr;
270 } 270 }
271 // clear PWA corrections table if any 271 // clear PWA corrections table if any
272 if (fPWACorrection) { 272 if (fPWACorrection) {
273 delete fPWACorrection; 273 delete fPWACorrection;
274 fPWACorrection = nullptr; 274 fPWACorrection = nullptr;
275 } 275 }
276 // create GS-table 276 // create GS-table
277 G4bool isElectron = true; 277 G4bool isElectron = true;
278 if (p->GetPDGCharge()>0.) { 278 if (p->GetPDGCharge()>0.) {
279 isElectron = false; 279 isElectron = false;
280 } 280 }
281 fGSTable = new G4GoudsmitSaundersonTable(i 281 fGSTable = new G4GoudsmitSaundersonTable(isElectron);
282 // G4GSTable will be initialised: 282 // G4GSTable will be initialised:
283 // - Screened-Rutherford DCS based GS angu 283 // - Screened-Rutherford DCS based GS angular distributions will be loaded only if they are not there yet
284 // - Mott-correction will be initialised i 284 // - Mott-correction will be initialised if Mott-correction was requested to be used
285 fGSTable->SetOptionMottCorrection(fIsUseMo 285 fGSTable->SetOptionMottCorrection(fIsUseMottCorrection);
286 // - set PWA correction (correction to int 286 // - set PWA correction (correction to integrated quantites from Dirac-PWA)
287 fGSTable->SetOptionPWACorrection(fIsUsePWA 287 fGSTable->SetOptionPWACorrection(fIsUsePWACorrection);
288 // init 288 // init
289 fGSTable->Initialise(LowEnergyLimit(),High 289 fGSTable->Initialise(LowEnergyLimit(),HighEnergyLimit());
290 // create PWA corrections table if it was 290 // create PWA corrections table if it was requested (and not disactivated because active Mott-correction)
291 if (fIsUsePWACorrection) { 291 if (fIsUsePWACorrection) {
292 fPWACorrection = new G4GSPWACorrections( 292 fPWACorrection = new G4GSPWACorrections(isElectron);
293 fPWACorrection->Initialise(); 293 fPWACorrection->Initialise();
294 } 294 }
295 } 295 }
296 fParticleChange = GetParticleChangeForMSC(p) 296 fParticleChange = GetParticleChangeForMSC(p);
297 } 297 }
298 298
299 299
300 void G4GoudsmitSaundersonMscModel::InitialiseL 300 void G4GoudsmitSaundersonMscModel::InitialiseLocal(const G4ParticleDefinition*, G4VEmModel* masterModel) {
301 fGSTable = static_cast<G4Goud 301 fGSTable = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetGSTable();
302 fIsUseMottCorrection = static_cast<G4Goud 302 fIsUseMottCorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetOptionMottCorrection();
303 fIsUsePWACorrection = static_cast<G4Goud 303 fIsUsePWACorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetOptionPWACorrection();
304 fPWACorrection = static_cast<G4Goud 304 fPWACorrection = static_cast<G4GoudsmitSaundersonMscModel*>(masterModel)->GetPWACorrection();
305 } 305 }
306 306
307 307
308 // computes macroscopic first transport cross 308 // computes macroscopic first transport cross section: used only in testing not during mc transport
309 G4double G4GoudsmitSaundersonMscModel::CrossSe 309 G4double G4GoudsmitSaundersonMscModel::CrossSectionPerVolume(const G4Material* mat,
310 const 310 const G4ParticleDefinition*,
311 G4dou 311 G4double kineticEnergy,
312 G4dou 312 G4double,
313 G4dou 313 G4double) {
314 G4double xsecTr1 = 0.; // cross section per 314 G4double xsecTr1 = 0.; // cross section per volume i.e. macroscopic 1st transport cross section
315 // 315 //
316 fLambda0 = 0.0; // elastic mean free path 316 fLambda0 = 0.0; // elastic mean free path
317 fLambda1 = 0.0; // first transport mean free 317 fLambda1 = 0.0; // first transport mean free path
318 fScrA = 0.0; // screening parameter 318 fScrA = 0.0; // screening parameter
319 fG1 = 0.0; // first transport coef. 319 fG1 = 0.0; // first transport coef.
320 // use Moliere's screening (with Mott-corret 320 // use Moliere's screening (with Mott-corretion if it was requested)
321 G4double efEnergy = std::max(kineticEnergy, 321 G4double efEnergy = std::max(kineticEnergy, 10.*CLHEP::eV);
322 // total mometum square 322 // total mometum square
323 G4double pt2 = efEnergy*(efEnergy+2.0*el 323 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
324 // beta square 324 // beta square
325 G4double beta2 = pt2/(pt2+electron_mass_c2 325 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
326 // current material index 326 // current material index
327 G4int matindx = (G4int)mat->GetIndex(); << 327 G4int matindx = mat->GetIndex();
328 // Moliere's b_c 328 // Moliere's b_c
329 G4double bc = fGSTable->GetMoliereBc(ma 329 G4double bc = fGSTable->GetMoliereBc(matindx);
330 // get the Mott-correcton factors if Mott-co 330 // get the Mott-correcton factors if Mott-correcton was requested by the user
331 fMCtoScrA = 1.0; 331 fMCtoScrA = 1.0;
332 fMCtoQ1 = 1.0; 332 fMCtoQ1 = 1.0;
333 fMCtoG2PerG1 = 1.0; 333 fMCtoG2PerG1 = 1.0;
334 G4double scpCor = 1.0; 334 G4double scpCor = 1.0;
335 if (fIsUseMottCorrection) { 335 if (fIsUseMottCorrection) {
336 fGSTable->GetMottCorrectionFactors(G4Log(e 336 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
337 // ! no scattering power correction since 337 // ! no scattering power correction since the current couple is not set before this interface method is called
338 // scpCor = fGSTable->ComputeScatteringPow 338 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
339 } else if (fIsUsePWACorrection) { 339 } else if (fIsUsePWACorrection) {
340 fPWACorrection->GetPWACorrectionFactors(G4 340 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
341 // scpCor = fGSTable->ComputeScatteringPow 341 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
342 } 342 }
343 // screening parameter: 343 // screening parameter:
344 // - if Mott-corretioncorrection: the Screen 344 // - if Mott-corretioncorrection: the Screened-Rutherford times Mott-corretion DCS with this
345 // screening parameter gives back the (els 345 // screening parameter gives back the (elsepa) PWA first transport cross section
346 // - if PWA correction: he Screened-Rutherfo 346 // - if PWA correction: he Screened-Rutherford DCS with this screening parameter
347 // gives back the (elsepa) PWA first trans 347 // gives back the (elsepa) PWA first transport cross section
348 fScrA = fGSTable->GetMoliereXc2(matindx)/ 348 fScrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*fMCtoScrA;
349 // elastic mean free path in Geant4 internal 349 // elastic mean free path in Geant4 internal lenght units: the neglected (1+screening parameter) term is corrected
350 // (if Mott-corretion: the corrected screeni 350 // (if Mott-corretion: the corrected screening parameter is used for this (1+A) correction + Moliere b_c is also
351 // corrected with the screening parameter co 351 // corrected with the screening parameter correction)
352 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scp 352 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scpCor;
353 // first transport coefficient (if Mott-corr 353 // first transport coefficient (if Mott-corretion: the corrected screening parameter is used (it will be fully
354 // consistent with the one used during the p 354 // consistent with the one used during the pre-computation of the Mott-correted GS angular distributions))
355 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/ 355 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0);
356 // first transport mean free path 356 // first transport mean free path
357 fLambda1 = fLambda0/fG1; 357 fLambda1 = fLambda0/fG1;
358 xsecTr1 = 1./fLambda1; 358 xsecTr1 = 1./fLambda1;
359 return xsecTr1; 359 return xsecTr1;
360 } 360 }
361 361
362 362
363 // gives back the first transport mean free pa 363 // gives back the first transport mean free path in internal G4 units
364 G4double 364 G4double
365 G4GoudsmitSaundersonMscModel::GetTransportMean 365 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(const G4ParticleDefinition* /*partdef*/,
366 366 G4double kineticEnergy) {
367 // kinetic energy is assumed to be in Geant4 367 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
368 G4double efEnergy = kineticEnergy; 368 G4double efEnergy = kineticEnergy;
369 // 369 //
370 const G4Material* mat = currentCouple->GetM 370 const G4Material* mat = currentCouple->GetMaterial();
371 // 371 //
372 fLambda0 = 0.0; // elastic mean free path 372 fLambda0 = 0.0; // elastic mean free path
373 fLambda1 = 0.0; // first transport mean free 373 fLambda1 = 0.0; // first transport mean free path
374 fScrA = 0.0; // screening parameter 374 fScrA = 0.0; // screening parameter
375 fG1 = 0.0; // first transport coef. 375 fG1 = 0.0; // first transport coef.
376 376
377 // use Moliere's screening (with Mott-corret 377 // use Moliere's screening (with Mott-corretion if it was requested)
378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* 378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*CLHEP::eV;
379 // total mometum square 379 // total mometum square
380 G4double pt2 = efEnergy*(efEnergy+2.0*el 380 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
381 // beta square 381 // beta square
382 G4double beta2 = pt2/(pt2+electron_mass_c2 382 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
383 // current material index 383 // current material index
384 G4int matindx = (G4int)mat->GetIndex(); << 384 G4int matindx = mat->GetIndex();
385 // Moliere's b_c 385 // Moliere's b_c
386 G4double bc = fGSTable->GetMoliereBc(ma 386 G4double bc = fGSTable->GetMoliereBc(matindx);
387 // get the Mott-correcton factors if Mott-co 387 // get the Mott-correcton factors if Mott-correcton was requested by the user
388 fMCtoScrA = 1.0; 388 fMCtoScrA = 1.0;
389 fMCtoQ1 = 1.0; 389 fMCtoQ1 = 1.0;
390 fMCtoG2PerG1 = 1.0; 390 fMCtoG2PerG1 = 1.0;
391 G4double scpCor = 1.0; 391 G4double scpCor = 1.0;
392 if (fIsUseMottCorrection) { 392 if (fIsUseMottCorrection) {
393 fGSTable->GetMottCorrectionFactors(G4Log(e 393 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
394 scpCor = fGSTable->ComputeScatteringPowerC 394 scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
395 } else if (fIsUsePWACorrection) { 395 } else if (fIsUsePWACorrection) {
396 fPWACorrection->GetPWACorrectionFactors(G4 396 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, fMCtoScrA, fMCtoQ1, fMCtoG2PerG1);
397 // scpCor = fGSTable->ComputeScatteringPow 397 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
398 } 398 }
399 // screening parameter: 399 // screening parameter:
400 // - if Mott-corretioncorrection: the Screen 400 // - if Mott-corretioncorrection: the Screened-Rutherford times Mott-corretion DCS with this
401 // screening parameter gives back the (els 401 // screening parameter gives back the (elsepa) PWA first transport cross section
402 // - if PWA correction: he Screened-Rutherfo 402 // - if PWA correction: he Screened-Rutherford DCS with this screening parameter
403 // gives back the (elsepa) PWA first trans 403 // gives back the (elsepa) PWA first transport cross section
404 fScrA = fGSTable->GetMoliereXc2(matindx)/ 404 fScrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*fMCtoScrA;
405 // elastic mean free path in Geant4 internal 405 // elastic mean free path in Geant4 internal lenght units: the neglected (1+screening parameter) term is corrected
406 // (if Mott-corretion: the corrected screeni 406 // (if Mott-corretion: the corrected screening parameter is used for this (1+A) correction + Moliere b_c is also
407 // corrected with the screening parameter co 407 // corrected with the screening parameter correction)
408 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scp 408 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scpCor;
409 // first transport coefficient (if Mott-corr 409 // first transport coefficient (if Mott-corretion: the corrected screening parameter is used (it will be fully
410 // consistent with the one used during the p 410 // consistent with the one used during the pre-computation of the Mott-correted GS angular distributions))
411 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/ 411 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0);
412 // first transport mean free path 412 // first transport mean free path
413 fLambda1 = fLambda0/fG1; 413 fLambda1 = fLambda0/fG1;
414 414
415 return fLambda1; 415 return fLambda1;
416 } 416 }
417 417
418 418
419 G4double 419 G4double
420 G4GoudsmitSaundersonMscModel::GetTransportMean 420 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly(const G4ParticleDefinition* /*partdef*/,
421 421 G4double kineticEnergy) {
422 // kinetic energy is assumed to be in Geant4 422 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
423 G4double efEnergy = kineticEnergy; 423 G4double efEnergy = kineticEnergy;
424 // 424 //
425 const G4Material* mat = currentCouple->GetM 425 const G4Material* mat = currentCouple->GetMaterial();
426 // 426 //
427 G4double lambda0 = 0.0; // elastc mean free 427 G4double lambda0 = 0.0; // elastc mean free path
428 G4double lambda1 = 0.0; // first transport m 428 G4double lambda1 = 0.0; // first transport mean free path
429 G4double scrA = 0.0; // screening paramet 429 G4double scrA = 0.0; // screening parametr
430 G4double g1 = 0.0; // first transport m 430 G4double g1 = 0.0; // first transport mean free path
431 431
432 // use Moliere's screening (with Mott-corret 432 // use Moliere's screening (with Mott-corretion if it was requested)
433 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* 433 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*CLHEP::eV;
434 // total mometum square in Geant4 internal e 434 // total mometum square in Geant4 internal energy2 units which is MeV2
435 G4double pt2 = efEnergy*(efEnergy+2.0*el 435 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
436 G4double beta2 = pt2/(pt2+electron_mass_c2 436 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
437 G4int matindx = (G4int)mat->GetIndex(); << 437 G4int matindx = mat->GetIndex();
438 G4double bc = fGSTable->GetMoliereBc(ma 438 G4double bc = fGSTable->GetMoliereBc(matindx);
439 // get the Mott-correcton factors if Mott-co 439 // get the Mott-correcton factors if Mott-correcton was requested by the user
440 G4double mctoScrA = 1.0; 440 G4double mctoScrA = 1.0;
441 G4double mctoQ1 = 1.0; 441 G4double mctoQ1 = 1.0;
442 G4double mctoG2PerG1 = 1.0; 442 G4double mctoG2PerG1 = 1.0;
443 G4double scpCor = 1.0; 443 G4double scpCor = 1.0;
444 if (fIsUseMottCorrection) { 444 if (fIsUseMottCorrection) {
445 fGSTable->GetMottCorrectionFactors(G4Log(e 445 fGSTable->GetMottCorrectionFactors(G4Log(efEnergy), beta2, matindx, mctoScrA, mctoQ1, mctoG2PerG1);
446 scpCor = fGSTable->ComputeScatteringPowerC 446 scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
447 } else if (fIsUsePWACorrection) { 447 } else if (fIsUsePWACorrection) {
448 fPWACorrection->GetPWACorrectionFactors(G4 448 fPWACorrection->GetPWACorrectionFactors(G4Log(efEnergy), beta2, matindx, mctoScrA, mctoQ1, mctoG2PerG1);
449 // scpCor = fGSTable->ComputeScatteringPow 449 // scpCor = fGSTable->ComputeScatteringPowerCorrection(currentCouple, efEnergy);
450 } 450 }
451 scrA = fGSTable->GetMoliereXc2(matindx)/( 451 scrA = fGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc)*mctoScrA;
452 // total elastic mean free path in Geant4 in 452 // total elastic mean free path in Geant4 internal lenght units
453 lambda0 = beta2*(1.+scrA)*mctoScrA/bc/scpCor 453 lambda0 = beta2*(1.+scrA)*mctoScrA/bc/scpCor;
454 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scr 454 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scrA+1.0)-1.0);
455 lambda1 = lambda0/g1; 455 lambda1 = lambda0/g1;
456 456
457 return lambda1; 457 return lambda1;
458 } 458 }
459 459
460 460
461 void G4GoudsmitSaundersonMscModel::StartTracki 461 void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track) {
462 SetParticle(track->GetDynamicParticle()->Get 462 SetParticle(track->GetDynamicParticle()->GetDefinition());
463 firstStep = true; 463 firstStep = true;
464 tlimit = tgeom = rangeinit = geombig; 464 tlimit = tgeom = rangeinit = geombig;
465 rangeinit = 1.e+21; 465 rangeinit = 1.e+21;
466 } 466 }
467 467
468 468
469 G4double G4GoudsmitSaundersonMscModel::Compute 469 G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
470 470 G4double& currentMinimalStep) {
471 G4double skindepth = 0.; 471 G4double skindepth = 0.;
472 // 472 //
473 const G4DynamicParticle* dp = track.GetDynam 473 const G4DynamicParticle* dp = track.GetDynamicParticle();
474 G4StepPoint* sp = track.GetStep( 474 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
475 G4StepStatus stepStatus = sp->GetStepSta 475 G4StepStatus stepStatus = sp->GetStepStatus();
476 currentCouple = track.GetMater 476 currentCouple = track.GetMaterialCutsCouple();
477 SetCurrentCouple(currentCouple); 477 SetCurrentCouple(currentCouple);
478 currentMaterialIndex = (G4int)current << 478 currentMaterialIndex = currentCouple->GetMaterial()->GetIndex();
479 currentKinEnergy = dp->GetKinetic 479 currentKinEnergy = dp->GetKineticEnergy();
480 currentRange = GetRange(parti 480 currentRange = GetRange(particle,currentKinEnergy,currentCouple,
481 dp->G 481 dp->GetLogKineticEnergy());
482 // elastic and first transport mfp, screenin 482 // elastic and first transport mfp, screening parameter and G1 are also set
483 // (Mott-correction will be used if it was r 483 // (Mott-correction will be used if it was requested by the user)
484 fLambda1 = GetTransportMeanFreePath(particle 484 fLambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
485 // Set initial values: 485 // Set initial values:
486 // : lengths are initialised to currentMini 486 // : lengths are initialised to currentMinimalStep which is the true, minimum
487 // step length from all other physics 487 // step length from all other physics
488 fTheTrueStepLenght = currentMinimalStep; 488 fTheTrueStepLenght = currentMinimalStep;
489 fTheTransportDistance = currentMinimalStep; 489 fTheTransportDistance = currentMinimalStep;
490 fTheZPathLenght = currentMinimalStep; 490 fTheZPathLenght = currentMinimalStep; // will need to be converted
491 fTheDisplacementVector.set(0.,0.,0.); 491 fTheDisplacementVector.set(0.,0.,0.);
492 fTheNewDirection.set(0.,0.,1.); 492 fTheNewDirection.set(0.,0.,1.);
493 493
494 // Can everything be done in the step limit 494 // Can everything be done in the step limit phase ?
495 fIsEverythingWasDone = false; 495 fIsEverythingWasDone = false;
496 // Multiple scattering needs to be sample ? 496 // Multiple scattering needs to be sample ?
497 fIsMultipleSacettring = false; 497 fIsMultipleSacettring = false;
498 // Single scattering needs to be sample ? 498 // Single scattering needs to be sample ?
499 fIsSingleScattering = false; 499 fIsSingleScattering = false;
500 // Was zero deflection in multiple scatterin 500 // Was zero deflection in multiple scattering sampling ?
501 fIsNoScatteringInMSC = false; 501 fIsNoScatteringInMSC = false;
502 // Do not care about displacement in MSC sam 502 // Do not care about displacement in MSC sampling
503 // ( used only in the case of gIsOptimizatio 503 // ( used only in the case of gIsOptimizationOn = true)
504 fIsNoDisplace = false; 504 fIsNoDisplace = false;
505 // get pre-step point safety 505 // get pre-step point safety
506 presafety = sp->GetSafety(); 506 presafety = sp->GetSafety();
507 // 507 //
508 fZeff = currentCouple->GetMaterial()->GetIon 508 fZeff = currentCouple->GetMaterial()->GetIonisation()->GetZeffective();
509 // distance will take into account max-fluct 509 // distance will take into account max-fluct.
510 G4double distance = currentRange; 510 G4double distance = currentRange;
511 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZe 511 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZeff));
512 // 512 //
513 // Possible optimization : if the distance i 513 // Possible optimization : if the distance is samller than the safety -> the
514 // particle will never leave this volume -> 514 // particle will never leave this volume -> dispalcement
515 // as the effect of multiple elastic scatter 515 // as the effect of multiple elastic scattering can be skipped
516 // Important : this optimization can cause p 516 // Important : this optimization can cause problems if one does scoring
517 // in a bigger volume since MSC won't be don 517 // in a bigger volume since MSC won't be done deep inside the volume when
518 // distance < safety so don't use optimized- 518 // distance < safety so don't use optimized-mode in such case.
519 if (gIsOptimizationOn && (distance<presafety 519 if (gIsOptimizationOn && (distance<presafety)) {
520 // Indicate that we need to do MSC after 520 // Indicate that we need to do MSC after transportation and no dispalcement.
521 fIsMultipleSacettring = true; 521 fIsMultipleSacettring = true;
522 fIsNoDisplace = true; 522 fIsNoDisplace = true;
523 } else if (steppingAlgorithm==fUseDistanceTo 523 } else if (steppingAlgorithm==fUseDistanceToBoundary) {
524 //Compute geomlimit (and presafety) : 524 //Compute geomlimit (and presafety) :
525 // - geomlimit will be: 525 // - geomlimit will be:
526 // == the straight line distance to the 526 // == the straight line distance to the boundary if currentRange is
527 // longer than that 527 // longer than that
528 // == a big value [geombig = 1.e50*mm] 528 // == a big value [geombig = 1.e50*mm] if currentRange is shorter than
529 // the straight line distance to the 529 // the straight line distance to the boundary
530 // - presafety will be updated as well 530 // - presafety will be updated as well
531 // So the particle can travell 'gemlimit' 531 // So the particle can travell 'gemlimit' distance (along a straight
532 // line!) in its current direction: 532 // line!) in its current direction:
533 // (1) before reaching a boundary (geomli 533 // (1) before reaching a boundary (geomlimit < geombig) OR
534 // (2) before reaching its current range 534 // (2) before reaching its current range (geomlimit == geombig)
535 geomlimit = ComputeGeomLimit(track, presaf 535 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
536 // Record that the particle is on a bounda 536 // Record that the particle is on a boundary
537 if ( (stepStatus==fGeomBoundary) || (stepS 537 if ( (stepStatus==fGeomBoundary) || (stepStatus==fUndefined && presafety==0.0)) {
538 fIsWasOnBoundary = true; 538 fIsWasOnBoundary = true;
539 } 539 }
540 // Set skin depth = skin x elastic_mean_fr 540 // Set skin depth = skin x elastic_mean_free_path
541 skindepth = skin*fLambda0; 541 skindepth = skin*fLambda0;
542 // Init the flag that indicates that the p 542 // Init the flag that indicates that the particle are within a skindepth
543 // distance from a boundary 543 // distance from a boundary
544 fIsInsideSkin = false; 544 fIsInsideSkin = false;
545 // Check if we can try Single Scattering b 545 // Check if we can try Single Scattering because we are within skindepth
546 // distance from/to a boundary OR the curr 546 // distance from/to a boundary OR the current minimum true-step-length is
547 // shorter than skindepth. NOTICE: the lat 547 // shorter than skindepth. NOTICE: the latest has only efficieny reasons
548 // because the MSC angular sampling is fin 548 // because the MSC angular sampling is fine for any short steps but much
549 // faster to try single scattering in case 549 // faster to try single scattering in case of short steps.
550 if ((stepStatus==fGeomBoundary) || (presaf 550 if ((stepStatus==fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) {
551 // check if we are within skindepth dist 551 // check if we are within skindepth distance from a boundary
552 if ((stepStatus == fGeomBoundary) || (pr 552 if ((stepStatus == fGeomBoundary) || (presafety < skindepth)) {
553 fIsInsideSkin = true; 553 fIsInsideSkin = true;
554 fIsWasOnBoundary = true; 554 fIsWasOnBoundary = true;
555 } 555 }
556 //Try single scattering: 556 //Try single scattering:
557 // - sample distance to next single scat 557 // - sample distance to next single scattering interaction (sslimit)
558 // - compare to current minimum length 558 // - compare to current minimum length
559 // == if sslimit is the shorter: 559 // == if sslimit is the shorter:
560 // - set the step length to ssl 560 // - set the step length to sslimit
561 // - indicate that single scatt 561 // - indicate that single scattering needs to be done
562 // == else : nothing to do 562 // == else : nothing to do
563 //- in both cases, the step length was v 563 //- in both cases, the step length was very short so geometrical and
564 // true path length are the same 564 // true path length are the same
565 G4double sslimit = -1.*fLambda0*G4Log(G4 565 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand());
566 // compare to current minimum step lengt 566 // compare to current minimum step length
567 if (sslimit<fTheTrueStepLenght) { 567 if (sslimit<fTheTrueStepLenght) {
568 fTheTrueStepLenght = sslimit; 568 fTheTrueStepLenght = sslimit;
569 fIsSingleScattering = true; 569 fIsSingleScattering = true;
570 } 570 }
571 // short step -> true step length equal 571 // short step -> true step length equal to geometrical path length
572 fTheZPathLenght = fTheTrueStepLeng 572 fTheZPathLenght = fTheTrueStepLenght;
573 // Set taht everything is done in step-l 573 // Set taht everything is done in step-limit phase so no MSC call
574 // We will check if we need to perform t 574 // We will check if we need to perform the single-scattering angular
575 // sampling i.e. if single elastic scatt 575 // sampling i.e. if single elastic scattering was the winer!
576 fIsEverythingWasDone = true; 576 fIsEverythingWasDone = true;
577 } else { 577 } else {
578 // After checking we know that we cannot 578 // After checking we know that we cannot try single scattering so we will
579 // need to make an MSC step 579 // need to make an MSC step
580 // Indicate that we need to make and MSC 580 // Indicate that we need to make and MSC step. We do not check if we can
581 // do it now i.e. if presafety>final_tru 581 // do it now i.e. if presafety>final_true_step_length so we let the
582 // fIsEverythingWasDone = false which in 582 // fIsEverythingWasDone = false which indicates that we will perform
583 // MSC after transportation. 583 // MSC after transportation.
584 fIsMultipleSacettring = true; 584 fIsMultipleSacettring = true;
585 // Init the first-real-step falg: it wil 585 // Init the first-real-step falg: it will indicate if we do the first
586 // non-single scattering step in this vo 586 // non-single scattering step in this volume with this particle
587 fIsFirstRealStep = false; 587 fIsFirstRealStep = false;
588 // If previously the partcile was on bou 588 // If previously the partcile was on boundary it was within skin as
589 // well. When it is not within skin anym 589 // well. When it is not within skin anymore it has just left the skin
590 // so we make the first real MSC step wi 590 // so we make the first real MSC step with the particle.
591 if (fIsWasOnBoundary && !fIsInsideSkin) 591 if (fIsWasOnBoundary && !fIsInsideSkin) {
592 // reset the 'was on boundary' indicat 592 // reset the 'was on boundary' indicator flag
593 fIsWasOnBoundary = false; 593 fIsWasOnBoundary = false;
594 fIsFirstRealStep = true; 594 fIsFirstRealStep = true;
595 } 595 }
596 // If this is the first-real msc step (t 596 // If this is the first-real msc step (the partcile has just left the
597 // skin) or this is the first step with 597 // skin) or this is the first step with the particle (was born or
598 // primary): 598 // primary):
599 // - set the initial range that will b 599 // - set the initial range that will be used later to limit its step
600 // (only in this volume, because aft 600 // (only in this volume, because after boundary crossing at the
601 // first-real MSC step we will reset 601 // first-real MSC step we will reset)
602 // - don't let the partcile to cross th 602 // - don't let the partcile to cross the volume just in one step
603 if (firstStep || fIsFirstRealStep || ran 603 if (firstStep || fIsFirstRealStep || rangeinit>1.e+20) {
604 rangeinit = currentRange; 604 rangeinit = currentRange;
605 // If geomlimit < geombig than the par 605 // If geomlimit < geombig than the particle might reach the boundary
606 // along its initial direction before 606 // along its initial direction before losing its energy (in this step)
607 // Otherwise we can be sure that the p 607 // Otherwise we can be sure that the particle will lose it energy
608 // before reaching the boundary along 608 // before reaching the boundary along a starigth line so there is no
609 // geometrical limit appalied. [Howeve 609 // geometrical limit appalied. [However, tgeom is set only in the
610 // first or the first-real MSC step. A 610 // first or the first-real MSC step. After the first or first real
611 // MSC step the direction will change 611 // MSC step the direction will change tgeom won't guaranty anything!
612 // But we will try to end up within sk 612 // But we will try to end up within skindepth from the boundary using
613 // the actual value of geomlimit(See l 613 // the actual value of geomlimit(See later at step reduction close to
614 // boundary).] 614 // boundary).]
615 if (geomlimit<geombig) { 615 if (geomlimit<geombig) {
616 // transfrom straight line distance 616 // transfrom straight line distance to the boundary to real step
617 // length based on the mean values ( 617 // length based on the mean values (using the prestep point
618 // first-transport mean free path i. 618 // first-transport mean free path i.e. no energy loss correction)
619 if ((1.-geomlimit/fLambda1)> 0.) { 619 if ((1.-geomlimit/fLambda1)> 0.) {
620 geomlimit = -fLambda1*G4Log(1.-geo 620 geomlimit = -fLambda1*G4Log(1.-geomlimit/fLambda1);
621 } 621 }
622 // the 2-different case that could l 622 // the 2-different case that could lead us here
623 if (firstStep) { 623 if (firstStep) {
624 tgeom = 2.*geomlimit/facgeom; 624 tgeom = 2.*geomlimit/facgeom;
625 } else { 625 } else {
626 tgeom = geomlimit/facgeom; 626 tgeom = geomlimit/facgeom;
627 } 627 }
628 } else { 628 } else {
629 tgeom = geombig; 629 tgeom = geombig;
630 } 630 }
631 } 631 }
632 // True step length limit from range fac 632 // True step length limit from range factor. Noteice, that the initial
633 // range is used that was set at the fir 633 // range is used that was set at the first step or first-real MSC step
634 // in this volume with this particle. 634 // in this volume with this particle.
635 tlimit = facrange*rangeinit; 635 tlimit = facrange*rangeinit;
636 // Take the minimum of the true step len 636 // Take the minimum of the true step length limits coming from
637 // geometrical constraint or range-facto 637 // geometrical constraint or range-factor limitation
638 tlimit = std::min(tlimit,tgeom); 638 tlimit = std::min(tlimit,tgeom);
639 // Step reduction close to boundary: we 639 // Step reduction close to boundary: we try to end up within skindepth
640 // from the boundary ( Notice: in case o 640 // from the boundary ( Notice: in case of mag. field it might not work
641 // because geomlimit is the straigth lin 641 // because geomlimit is the straigth line distance to the boundary in
642 // the currect direction (if geomlimit<g 642 // the currect direction (if geomlimit<geombig) and mag. field can
643 // change the initial direction. So te p 643 // change the initial direction. So te particle might hit some boundary
644 // before in a different direction. Howe 644 // before in a different direction. However, here we restrict the true
645 // path length to this (straight line) l 645 // path length to this (straight line) lenght so the corresponding
646 // transport distance (straight line) wi 646 // transport distance (straight line) will be even shorter than
647 // geomlimit-0.999*skindepth after the c 647 // geomlimit-0.999*skindepth after the change of true->geom.
648 if (geomlimit<geombig) { 648 if (geomlimit<geombig) {
649 tlimit = std::min(tlimit, geomlimit-0. 649 tlimit = std::min(tlimit, geomlimit-0.999*skindepth);
650 } 650 }
651 // randomize 1st step or 1st 'normal' st 651 // randomize 1st step or 1st 'normal' step in volume
652 if (firstStep || fIsFirstRealStep) { 652 if (firstStep || fIsFirstRealStep) {
653 fTheTrueStepLenght = std::min(fTheTrue 653 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
654 } else { 654 } else {
655 fTheTrueStepLenght = std::min(fTheTrue 655 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
656 } 656 }
657 } 657 }
658 } else if (steppingAlgorithm==fUseSafetyPlus 658 } else if (steppingAlgorithm==fUseSafetyPlus) { // THE ERROR_FREE stepping alg.
659 presafety = ComputeSafety(sp->GetPosition 659 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
660 geomlimit = presafety; 660 geomlimit = presafety;
661 // Set skin depth = skin x elastic_mean_fr 661 // Set skin depth = skin x elastic_mean_free_path
662 skindepth = skin*fLambda0; 662 skindepth = skin*fLambda0;
663 // Check if we can try Single Scattering b 663 // Check if we can try Single Scattering because we are within skindepth
664 // distance from/to a boundary OR the curr 664 // distance from/to a boundary OR the current minimum true-step-length is
665 // shorter than skindepth. NOTICE: the lat 665 // shorter than skindepth. NOTICE: the latest has only efficieny reasons
666 // because the MSC angular sampling is fin 666 // because the MSC angular sampling is fine for any short steps but much
667 // faster to try single scattering in case 667 // faster to try single scattering in case of short steps.
668 if ((stepStatus==fGeomBoundary) || (presaf 668 if ((stepStatus==fGeomBoundary) || (presafety<skindepth) || (fTheTrueStepLenght<skindepth)) {
669 //Try single scattering: 669 //Try single scattering:
670 // - sample distance to next single scat 670 // - sample distance to next single scattering interaction (sslimit)
671 // - compare to current minimum length 671 // - compare to current minimum length
672 // == if sslimit is the shorter: 672 // == if sslimit is the shorter:
673 // - set the step length to ssl 673 // - set the step length to sslimit
674 // - indicate that single scatt 674 // - indicate that single scattering needs to be done
675 // == else : nothing to do 675 // == else : nothing to do
676 //- in both cases, the step length was v 676 //- in both cases, the step length was very short so geometrical and
677 // true path length are the same 677 // true path length are the same
678 G4double sslimit = -1.*fLambda0*G4Log(G4 678 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand());
679 // compare to current minimum step lengt 679 // compare to current minimum step length
680 if (sslimit<fTheTrueStepLenght) { 680 if (sslimit<fTheTrueStepLenght) {
681 fTheTrueStepLenght = sslimit; 681 fTheTrueStepLenght = sslimit;
682 fIsSingleScattering = true; 682 fIsSingleScattering = true;
683 } 683 }
684 // short step -> true step length equal 684 // short step -> true step length equal to geometrical path length
685 fTheZPathLenght = fTheTrueStepLeng 685 fTheZPathLenght = fTheTrueStepLenght;
686 // Set taht everything is done in step-l 686 // Set taht everything is done in step-limit phase so no MSC call
687 // We will check if we need to perform t 687 // We will check if we need to perform the single-scattering angular
688 // sampling i.e. if single elastic scatt 688 // sampling i.e. if single elastic scattering was the winer!
689 fIsEverythingWasDone = true; 689 fIsEverythingWasDone = true;
690 } else { 690 } else {
691 // After checking we know that we cannot 691 // After checking we know that we cannot try single scattering so we will
692 // need to make an MSC step 692 // need to make an MSC step
693 // Indicate that we need to make and MSC 693 // Indicate that we need to make and MSC step.
694 fIsMultipleSacettring = true; 694 fIsMultipleSacettring = true;
695 fIsEverythingWasDone = true; 695 fIsEverythingWasDone = true;
696 // limit from range factor 696 // limit from range factor
697 fTheTrueStepLenght = std::min(fTheTru 697 fTheTrueStepLenght = std::min(fTheTrueStepLenght, facrange*currentRange);
698 // never let the particle go further tha 698 // never let the particle go further than the safety if we are out of the skin
699 // if we are here we are out of the skin 699 // if we are here we are out of the skin, presafety > 0.
700 if (fTheTrueStepLenght>presafety) { 700 if (fTheTrueStepLenght>presafety) {
701 fTheTrueStepLenght = std::min(fTheTrue 701 fTheTrueStepLenght = std::min(fTheTrueStepLenght, presafety);
702 } 702 }
703 // make sure that we are still within th 703 // make sure that we are still within the aplicability of condensed histry model
704 // i.e. true step length is not longer t 704 // i.e. true step length is not longer than first transport mean free path.
705 // We schould take into account energy l 705 // We schould take into account energy loss along 0.5x lambda_transport1
706 // step length as well. So let it 0.5 x 706 // step length as well. So let it 0.5 x lambda_transport1
707 fTheTrueStepLenght = std::min(fTheTrueSt 707 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*0.5);
708 } 708 }
709 } else { 709 } else {
710 // This is the default stepping algorithm: 710 // This is the default stepping algorithm: the fastest but the least
711 // accurate that corresponds to fUseSafety 711 // accurate that corresponds to fUseSafety in Urban model. Note, that GS
712 // model can handle any short steps so we 712 // model can handle any short steps so we do not need the minimum limits
713 // 713 //
714 // NO single scattering in case of skin or 714 // NO single scattering in case of skin or short steps (by defult the MSC
715 // model will be single or even no scatter 715 // model will be single or even no scattering in case of short steps
716 // compared to the elastic mean free path. 716 // compared to the elastic mean free path.)
717 // 717 //
718 // indicate that MSC needs to be done (alw 718 // indicate that MSC needs to be done (always and always after transportation)
719 fIsMultipleSacettring = true; 719 fIsMultipleSacettring = true;
720 if (stepStatus!=fGeomBoundary) { 720 if (stepStatus!=fGeomBoundary) {
721 presafety = ComputeSafety(sp->GetPositio 721 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
722 } 722 }
723 // Far from boundary-> in optimized mode d 723 // Far from boundary-> in optimized mode do not sample dispalcement.
724 if ((distance<presafety) && (gIsOptimizati 724 if ((distance<presafety) && (gIsOptimizationOn)) {
725 fIsNoDisplace = true; 725 fIsNoDisplace = true;
726 } else { 726 } else {
727 // Urban like 727 // Urban like
728 if (firstStep || (stepStatus==fGeomBound 728 if (firstStep || (stepStatus==fGeomBoundary) || rangeinit>1.e+20) {
729 rangeinit = currentRange; 729 rangeinit = currentRange;
730 fr = facrange; 730 fr = facrange;
731 // We don't use this: we won't converge to the 731 // We don't use this: we won't converge to the single scattering results with
732 // decreasing range-factor. 732 // decreasing range-factor.
733 // rangeinit = std::max(rangeinit 733 // rangeinit = std::max(rangeinit, fLambda1);
734 // if(fLambda1 > lambdalimit) { 734 // if(fLambda1 > lambdalimit) {
735 // fr *= (0.75+0.25*fLambda1/la 735 // fr *= (0.75+0.25*fLambda1/lambdalimit);
736 // } 736 // }
737 737
738 } 738 }
739 //step limit 739 //step limit
740 tlimit = std::max(fr*rangeinit, facsafet 740 tlimit = std::max(fr*rangeinit, facsafety*presafety);
741 // first step randomization 741 // first step randomization
742 if (firstStep || stepStatus==fGeomBounda 742 if (firstStep || stepStatus==fGeomBoundary) {
743 fTheTrueStepLenght = std::min(fTheTrue 743 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
744 } else { 744 } else {
745 fTheTrueStepLenght = std::min(fTheTrue 745 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
746 } 746 }
747 } 747 }
748 } 748 }
749 // 749 //
750 // unset first-step 750 // unset first-step
751 firstStep =false; 751 firstStep =false;
752 // performe single scattering, multiple scat 752 // performe single scattering, multiple scattering if this later can be done safely here
753 if (fIsEverythingWasDone) { 753 if (fIsEverythingWasDone) {
754 if (fIsSingleScattering) { 754 if (fIsSingleScattering) {
755 // sample single scattering 755 // sample single scattering
756 //G4double ekin = 0.5*(currentKinEnerg 756 //G4double ekin = 0.5*(currentKinEnergy + GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple));
757 G4double lekin = G4Log(currentKinEnergy 757 G4double lekin = G4Log(currentKinEnergy);
758 G4double pt2 = currentKinEnergy*(curr 758 G4double pt2 = currentKinEnergy*(currentKinEnergy+2.0*CLHEP::electron_mass_c2);
759 G4double beta2 = pt2/(pt2+CLHEP::electr 759 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
760 G4double cost = fGSTable->SingleScatte 760 G4double cost = fGSTable->SingleScattering(1., fScrA, lekin, beta2, currentMaterialIndex);
761 // protection 761 // protection
762 if (cost<-1.) cost = -1.; 762 if (cost<-1.) cost = -1.;
763 if (cost> 1.) cost = 1.; 763 if (cost> 1.) cost = 1.;
764 // compute sint 764 // compute sint
765 G4double dum = 1.-cost; 765 G4double dum = 1.-cost;
766 G4double sint = std::sqrt(dum*(2.-dum) 766 G4double sint = std::sqrt(dum*(2.-dum));
767 G4double phi = CLHEP::twopi*G4Uniform 767 G4double phi = CLHEP::twopi*G4UniformRand();
768 G4double sinPhi = std::sin(phi); 768 G4double sinPhi = std::sin(phi);
769 G4double cosPhi = std::cos(phi); 769 G4double cosPhi = std::cos(phi);
770 fTheNewDirection.set(sint*cosPhi,sint*si 770 fTheNewDirection.set(sint*cosPhi,sint*sinPhi,cost);
771 } else if (fIsMultipleSacettring) { 771 } else if (fIsMultipleSacettring) {
772 // sample multiple scattering 772 // sample multiple scattering
773 SampleMSC(); // fTheZPathLenght, fTheDis 773 SampleMSC(); // fTheZPathLenght, fTheDisplacementVector and fTheNewDirection will be set
774 } // and if single scattering but it was l 774 } // and if single scattering but it was longer => nothing to do
775 } //else { do nothing here but after transpo 775 } //else { do nothing here but after transportation
776 // 776 //
777 return ConvertTrueToGeom(fTheTrueStepLenght, 777 return ConvertTrueToGeom(fTheTrueStepLenght,currentMinimalStep);
778 } 778 }
779 779
780 780
781 G4double G4GoudsmitSaundersonMscModel::Compute 781 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double) {
782 // convert true ->geom 782 // convert true ->geom
783 // It is called from the step limitation Com 783 // It is called from the step limitation ComputeTruePathLengthLimit if
784 // !fIsEverythingWasDone but protect: 784 // !fIsEverythingWasDone but protect:
785 par1 = -1.; 785 par1 = -1.;
786 par2 = par3 = 0.; 786 par2 = par3 = 0.;
787 // if fIsEverythingWasDone = TRUE => fTheZP 787 // if fIsEverythingWasDone = TRUE => fTheZPathLenght is already set
788 // so return with the already known value 788 // so return with the already known value
789 // Otherwise: 789 // Otherwise:
790 if (!fIsEverythingWasDone) { 790 if (!fIsEverythingWasDone) {
791 // this correction needed to run MSC with 791 // this correction needed to run MSC with eIoni and eBrem inactivated
792 // and makes no harm for a normal run 792 // and makes no harm for a normal run
793 fTheTrueStepLenght = std::min(fTheTrueStep 793 fTheTrueStepLenght = std::min(fTheTrueStepLenght,currentRange);
794 // do the true -> geom transformation 794 // do the true -> geom transformation
795 fTheZPathLenght = fTheTrueStepLenght; 795 fTheZPathLenght = fTheTrueStepLenght;
796 // z = t for very small true-path-length 796 // z = t for very small true-path-length
797 if (fTheTrueStepLenght<tlimitminfix2) { 797 if (fTheTrueStepLenght<tlimitminfix2) {
798 return fTheZPathLenght; 798 return fTheZPathLenght;
799 } 799 }
800 G4double tau = fTheTrueStepLenght/fLambda1 800 G4double tau = fTheTrueStepLenght/fLambda1;
801 if (tau<=tausmall) { 801 if (tau<=tausmall) {
802 fTheZPathLenght = std::min(fTheTrueStepL 802 fTheZPathLenght = std::min(fTheTrueStepLenght, fLambda1);
803 } else if (fTheTrueStepLenght<currentRang 803 } else if (fTheTrueStepLenght<currentRange*dtrl) {
804 if (tau<taulim) fTheZPathLenght = fTheTr 804 if (tau<taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ;
805 else fTheZPathLenght = fLambd 805 else fTheZPathLenght = fLambda1*(1.-G4Exp(-tau));
806 } else if (currentKinEnergy<mass || fTheTr 806 } else if (currentKinEnergy<mass || fTheTrueStepLenght==currentRange) {
807 par1 = 1./currentRange ; // alpha =1 807 par1 = 1./currentRange ; // alpha =1/range_init for Ekin<mass
808 par2 = 1./(par1*fLambda1) ; // 1/(alpha 808 par2 = 1./(par1*fLambda1) ; // 1/(alphaxlambda01)
809 par3 = 1.+par2 ; // 1+1/ 809 par3 = 1.+par2 ; // 1+1/
810 if (fTheTrueStepLenght<currentRange) { 810 if (fTheTrueStepLenght<currentRange) {
811 fTheZPathLenght = 1./(par1*par3) * (1. 811 fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3));
812 } else { 812 } else {
813 fTheZPathLenght = 1./(par1*par3); 813 fTheZPathLenght = 1./(par1*par3);
814 } 814 }
815 } else { 815 } else {
816 G4double rfin = std::max(currentRange 816 G4double rfin = std::max(currentRange-fTheTrueStepLenght, 0.01*currentRange);
817 G4double T1 = GetEnergy(particle,rf 817 G4double T1 = GetEnergy(particle,rfin,currentCouple);
818 G4double lambda1 = GetTransportMeanFreeP 818 G4double lambda1 = GetTransportMeanFreePathOnly(particle,T1);
819 // 819 //
820 par1 = (fLambda1-lambda1)/(fLambda1*fThe 820 par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght); // alpha
821 par2 = 1./(par1*fLambda1); 821 par2 = 1./(par1*fLambda1);
822 par3 = 1.+par2 ; 822 par3 = 1.+par2 ;
823 G4Pow *g4calc = G4Pow::GetInstance(); 823 G4Pow *g4calc = G4Pow::GetInstance();
824 fTheZPathLenght = 1./(par1*par3) * (1.-g 824 fTheZPathLenght = 1./(par1*par3) * (1.-g4calc->powA(1.-par1*fTheTrueStepLenght,par3));
825 } 825 }
826 } 826 }
827 fTheZPathLenght = std::min(fTheZPathLenght, 827 fTheZPathLenght = std::min(fTheZPathLenght, fLambda1);
828 // 828 //
829 return fTheZPathLenght; 829 return fTheZPathLenght;
830 } 830 }
831 831
832 832
833 G4double G4GoudsmitSaundersonMscModel::Compute 833 G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength) {
834 // init 834 // init
835 fIsEndedUpOnBoundary = false; 835 fIsEndedUpOnBoundary = false;
836 // step defined other than transportation 836 // step defined other than transportation
837 if (geomStepLength==fTheZPathLenght) { 837 if (geomStepLength==fTheZPathLenght) {
838 return fTheTrueStepLenght; 838 return fTheTrueStepLenght;
839 } 839 }
840 // else :: 840 // else ::
841 // - set the flag that transportation was th 841 // - set the flag that transportation was the winer so DoNothin in DOIT !!
842 // - convert geom -> true by using the mean 842 // - convert geom -> true by using the mean value
843 fIsEndedUpOnBoundary = true; // OR LAST STEP 843 fIsEndedUpOnBoundary = true; // OR LAST STEP
844 fTheZPathLenght = geomStepLength; 844 fTheZPathLenght = geomStepLength;
845 // was a short single scattering step 845 // was a short single scattering step
846 if (fIsEverythingWasDone && !fIsMultipleSace 846 if (fIsEverythingWasDone && !fIsMultipleSacettring) {
847 fTheTrueStepLenght = geomStepLength; 847 fTheTrueStepLenght = geomStepLength;
848 return fTheTrueStepLenght; 848 return fTheTrueStepLenght;
849 } 849 }
850 // t = z for very small step 850 // t = z for very small step
851 if (geomStepLength<tlimitminfix2) { 851 if (geomStepLength<tlimitminfix2) {
852 fTheTrueStepLenght = geomStepLength; 852 fTheTrueStepLenght = geomStepLength;
853 // recalculation 853 // recalculation
854 } else { 854 } else {
855 G4double tlength = geomStepLength; 855 G4double tlength = geomStepLength;
856 if (geomStepLength>fLambda1*tausmall) { 856 if (geomStepLength>fLambda1*tausmall) {
857 if (par1< 0.) { 857 if (par1< 0.) {
858 tlength = -fLambda1*G4Log(1.-geomStepL 858 tlength = -fLambda1*G4Log(1.-geomStepLength/fLambda1) ;
859 } else { 859 } else {
860 if (par1*par3*geomStepLength<1.) { 860 if (par1*par3*geomStepLength<1.) {
861 G4Pow *g4calc = G4Pow::GetInstance() 861 G4Pow *g4calc = G4Pow::GetInstance();
862 tlength = (1.-g4calc->powA( 1.-par1* 862 tlength = (1.-g4calc->powA( 1.-par1*par3*geomStepLength,1./par3))/par1;
863 } else { 863 } else {
864 tlength = currentRange; 864 tlength = currentRange;
865 } 865 }
866 } 866 }
867 if (tlength<geomStepLength || tlength>fT 867 if (tlength<geomStepLength || tlength>fTheTrueStepLenght) {
868 tlength = geomStepLength; 868 tlength = geomStepLength;
869 } 869 }
870 } 870 }
871 fTheTrueStepLenght = tlength; 871 fTheTrueStepLenght = tlength;
872 } 872 }
873 // 873 //
874 return fTheTrueStepLenght; 874 return fTheTrueStepLenght;
875 } 875 }
876 876
877 G4ThreeVector& 877 G4ThreeVector&
878 G4GoudsmitSaundersonMscModel::SampleScattering 878 G4GoudsmitSaundersonMscModel::SampleScattering(const G4ThreeVector& oldDirection, G4double) {
879 if (steppingAlgorithm==fUseDistanceToBoundar 879 if (steppingAlgorithm==fUseDistanceToBoundary && fIsEverythingWasDone && fIsSingleScattering) {
880 // single scattering was and scattering ha 880 // single scattering was and scattering happend
881 fTheNewDirection.rotateUz(oldDirection); 881 fTheNewDirection.rotateUz(oldDirection);
882 fParticleChange->ProposeMomentumDirection( 882 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
883 return fTheDisplacementVector; 883 return fTheDisplacementVector;
884 } else if (steppingAlgorithm==fUseSafetyPlus 884 } else if (steppingAlgorithm==fUseSafetyPlus) { // error-free stepping
885 if (fIsEndedUpOnBoundary) { // do nothing 885 if (fIsEndedUpOnBoundary) { // do nothing on the boundary
886 return fTheDisplacementVector; 886 return fTheDisplacementVector;
887 } else if (fIsEverythingWasDone) { // evry 887 } else if (fIsEverythingWasDone) { // evrything is done if not optimizations case !!!
888 // check single scattering and see if it 888 // check single scattering and see if it happened
889 if (fIsSingleScattering) { 889 if (fIsSingleScattering) {
890 fTheNewDirection.rotateUz(oldDirection 890 fTheNewDirection.rotateUz(oldDirection);
891 fParticleChange->ProposeMomentumDirect 891 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
892 return fTheDisplacementVector; 892 return fTheDisplacementVector;
893 } 893 }
894 // check if multiple scattering happened 894 // check if multiple scattering happened and do things only if scattering was really happening
895 if (fIsMultipleSacettring && !fIsNoScatt 895 if (fIsMultipleSacettring && !fIsNoScatteringInMSC) {
896 fTheNewDirection.rotateUz(oldDirect 896 fTheNewDirection.rotateUz(oldDirection);
897 fTheDisplacementVector.rotateUz(old 897 fTheDisplacementVector.rotateUz(oldDirection);
898 fParticleChange->ProposeMomentumDir 898 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
899 } 899 }
900 // The only thing that could happen if w 900 // The only thing that could happen if we are here (fUseSafety and fIsEverythingWasDone)
901 // is that single scattering was tried 901 // is that single scattering was tried but did not win so scattering did not happen.
902 // So no displacement and no scattering 902 // So no displacement and no scattering
903 return fTheDisplacementVector; 903 return fTheDisplacementVector;
904 } 904 }
905 // 905 //
906 // The only thing that could still happen 906 // The only thing that could still happen with fUseSafetyPlus is that we are in the
907 // optimization branch: so sample MSC angl 907 // optimization branch: so sample MSC angle here (no displacement)
908 } 908 }
909 //else MSC needs to be done here 909 //else MSC needs to be done here
910 SampleMSC(); 910 SampleMSC();
911 if (!fIsNoScatteringInMSC) { 911 if (!fIsNoScatteringInMSC) {
912 fTheNewDirection.rotateUz(oldDirection); 912 fTheNewDirection.rotateUz(oldDirection);
913 fParticleChange->ProposeMomentumDirection( 913 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
914 if (!fIsNoDisplace) { 914 if (!fIsNoDisplace) {
915 fTheDisplacementVector.rotateUz(oldDirec 915 fTheDisplacementVector.rotateUz(oldDirection);
916 } 916 }
917 } 917 }
918 // 918 //
919 return fTheDisplacementVector; 919 return fTheDisplacementVector;
920 } 920 }
921 921
922 922
923 void G4GoudsmitSaundersonMscModel::SampleMSC() 923 void G4GoudsmitSaundersonMscModel::SampleMSC() {
924 fIsNoScatteringInMSC = false; 924 fIsNoScatteringInMSC = false;
925 // kinetic energy is assumed to be in Geant4 925 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
926 G4double kineticEnergy = currentKinEnergy; 926 G4double kineticEnergy = currentKinEnergy;
927 // 927 //
928 // Energy loss correction: 2 version 928 // Energy loss correction: 2 version
929 G4double eloss = 0.0; 929 G4double eloss = 0.0;
930 // if (fTheTrueStepLenght > currentRange*dtrl 930 // if (fTheTrueStepLenght > currentRange*dtrl) {
931 eloss = kineticEnergy - GetEnergy(particle,c 931 eloss = kineticEnergy - GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple);
932 // } else { 932 // } else {
933 // eloss = fTheTrueStepLenght*GetDEDX(parti 933 // eloss = fTheTrueStepLenght*GetDEDX(particle,kineticEnergy,currentCouple);
934 // } 934 // }
935 935
936 G4double tau = 0.;// = kineticEnergy/ele 936 G4double tau = 0.;// = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
937 G4double tau2 = 0.;// = tau*tau; 937 G4double tau2 = 0.;// = tau*tau;
938 G4double eps0 = 0.;// = eloss/kineticEnerg 938 G4double eps0 = 0.;// = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
939 G4double epsm = 0.;// = eloss/kineticEnerg 939 G4double epsm = 0.;// = eloss/kineticEnergy; // energy loss fraction to the mean step energy
940 940
941 // - init. 941 // - init.
942 G4double efEnergy = kineticEnergy; 942 G4double efEnergy = kineticEnergy;
943 G4double efStep = fTheTrueStepLenght; 943 G4double efStep = fTheTrueStepLenght;
944 944
945 G4double kineticEnergy0 = kineticEnergy; 945 G4double kineticEnergy0 = kineticEnergy;
946 if (gIsUseAccurate) { // - use accurate e 946 if (gIsUseAccurate) { // - use accurate energy loss correction
947 kineticEnergy -= 0.5*eloss; // mean ener 947 kineticEnergy -= 0.5*eloss; // mean energy along the full step
948 // other parameters for energy loss correc 948 // other parameters for energy loss corrections
949 tau = kineticEnergy/electron_m 949 tau = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
950 tau2 = tau*tau; 950 tau2 = tau*tau;
951 eps0 = eloss/kineticEnergy0; // 951 eps0 = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
952 epsm = eloss/kineticEnergy; // 952 epsm = eloss/kineticEnergy; // energy loss fraction to the mean step energy
953 953
954 efEnergy = kineticEnergy * (1.-epsm 954 efEnergy = kineticEnergy * (1.-epsm*epsm*(6.+10.*tau+5.*tau2)/(24.*tau2+48.*tau+72.));
955 G4double dum = 0.166666*(4.+tau*(6.+tau 955 G4double dum = 0.166666*(4.+tau*(6.+tau*(7.+tau*(4.+tau))))*(epsm/((tau+1.)*(tau+2.)))*(epsm/((tau+1.)*(tau+2.)));
956 efStep = fTheTrueStepLenght*(1.-d 956 efStep = fTheTrueStepLenght*(1.-dum);
957 } else { // - t 957 } else { // - take only mean energy
958 kineticEnergy -= 0.5*eloss; // mean ener 958 kineticEnergy -= 0.5*eloss; // mean energy along the full step
959 efEnergy = kineticEnergy; 959 efEnergy = kineticEnergy;
960 G4double factor = 1./(1.+0.9784671*kinetic 960 G4double factor = 1./(1.+0.9784671*kineticEnergy); //0.9784671 = 1/(2*m_e)
961 eps0 = eloss/kineticEnergy0; 961 eps0 = eloss/kineticEnergy0;
962 epsm = eps0/(1.-0.5*eps0); 962 epsm = eps0/(1.-0.5*eps0);
963 G4double temp = 0.3*(1 -factor*(1.-0.333 963 G4double temp = 0.3*(1 -factor*(1.-0.333333*factor))*eps0*eps0;
964 efStep = fTheTrueStepLenght*(1.+t 964 efStep = fTheTrueStepLenght*(1.+temp);
965 } 965 }
966 // 966 //
967 // compute elastic mfp, first transport mfp, 967 // compute elastic mfp, first transport mfp, screening parameter, and G1 (with Mott-correction
968 // if it was requested by the user) 968 // if it was requested by the user)
969 fLambda1 = GetTransportMeanFreePath(particle 969 fLambda1 = GetTransportMeanFreePath(particle, efEnergy);
970 // s/lambda_el 970 // s/lambda_el
971 G4double lambdan=0.; 971 G4double lambdan=0.;
972 if (fLambda0>0.0) { 972 if (fLambda0>0.0) {
973 lambdan=efStep/fLambda0; 973 lambdan=efStep/fLambda0;
974 } 974 }
975 if (lambdan<=1.0e-12) { 975 if (lambdan<=1.0e-12) {
976 if (fIsEverythingWasDone) { 976 if (fIsEverythingWasDone) {
977 fTheZPathLenght = fTheTrueStepLenght; 977 fTheZPathLenght = fTheTrueStepLenght;
978 } 978 }
979 fIsNoScatteringInMSC = true; 979 fIsNoScatteringInMSC = true;
980 return; 980 return;
981 } 981 }
982 // first moment: 2.* lambdan *scrA*((1.+scrA 982 // first moment: 2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
983 G4double Qn1 = lambdan *fG1; 983 G4double Qn1 = lambdan *fG1;
984 // sample scattering angles 984 // sample scattering angles
985 // new direction, relative to the orriginal 985 // new direction, relative to the orriginal one is in {uss,vss,wss}
986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, c 986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, cosTheta2 = 1.0, sinTheta2 = 0.0;
987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, c 987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, cosPhi2 = 1.0, sinPhi2 = 0.0;
988 G4double uss = 0.0, vss = 0.0, w 988 G4double uss = 0.0, vss = 0.0, wss = 1.0;
989 G4double x_coord = 0.0, y_coord = 0.0, z 989 G4double x_coord = 0.0, y_coord = 0.0, z_coord = 1.0;
990 G4double u2 = 0.0, v2 = 0.0; 990 G4double u2 = 0.0, v2 = 0.0;
991 // if we are above the upper grid limit with 991 // if we are above the upper grid limit with lambdaxG1=true-length/first-trans-mfp
992 // => izotropic distribution: lambG1_max =7. 992 // => izotropic distribution: lambG1_max =7.992 but set it to 7
993 if (0.5*Qn1 > 7.0){ 993 if (0.5*Qn1 > 7.0){
994 cosTheta1 = 1.-2.*G4UniformRand(); 994 cosTheta1 = 1.-2.*G4UniformRand();
995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+c 995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1));
996 cosTheta2 = 1.-2.*G4UniformRand(); 996 cosTheta2 = 1.-2.*G4UniformRand();
997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+c 997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+cosTheta2));
998 } else { 998 } else {
999 // sample 2 scattering cost1, sint1, cost 999 // sample 2 scattering cost1, sint1, cost2 and sint2 for half path
1000 G4double lekin = G4Log(efEnergy); 1000 G4double lekin = G4Log(efEnergy);
1001 G4double pt2 = efEnergy*(efEnergy+2.0 1001 G4double pt2 = efEnergy*(efEnergy+2.0*CLHEP::electron_mass_c2);
1002 G4double beta2 = pt2/(pt2+CLHEP::electr 1002 G4double beta2 = pt2/(pt2+CLHEP::electron_mass_c2*CLHEP::electron_mass_c2);
1003 // backup GS angular dtr pointer (kineti 1003 // backup GS angular dtr pointer (kinetic energy and delta index in case of Mott-correction)
1004 // if the first was an msc sampling (the 1004 // if the first was an msc sampling (the same will be used if the second is also an msc step)
1005 G4GoudsmitSaundersonTable::GSMSCAngularD 1005 G4GoudsmitSaundersonTable::GSMSCAngularDtr *gsDtr = nullptr;
1006 G4int mcEkinIdx = -1; 1006 G4int mcEkinIdx = -1;
1007 G4int mcDeltIdx = -1; 1007 G4int mcDeltIdx = -1;
1008 G4double transfPar = 0.; 1008 G4double transfPar = 0.;
1009 G4bool isMsc = fGSTable->Sampling(0.5*la 1009 G4bool isMsc = fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta1, sinTheta1, lekin, beta2,
1010 curren 1010 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar,
1011 true); 1011 true);
1012 fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, 1012 fGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta2, sinTheta2, lekin, beta2,
1013 currentMaterialIndex, 1013 currentMaterialIndex, &gsDtr, mcEkinIdx, mcDeltIdx, transfPar, !isMsc);
1014 if (cosTheta1+cosTheta2==2.) { // no sca 1014 if (cosTheta1+cosTheta2==2.) { // no scattering happened
1015 if (fIsEverythingWasDone) 1015 if (fIsEverythingWasDone)
1016 fTheZPathLenght = fTheTrueStepLeng 1016 fTheZPathLenght = fTheTrueStepLenght;
1017 fIsNoScatteringInMSC = true; 1017 fIsNoScatteringInMSC = true;
1018 return; 1018 return;
1019 } 1019 }
1020 } 1020 }
1021 // sample 2 azimuthal angles 1021 // sample 2 azimuthal angles
1022 G4double phi1 = CLHEP::twopi*G4UniformRand( 1022 G4double phi1 = CLHEP::twopi*G4UniformRand();
1023 sinPhi1 = std::sin(phi1); 1023 sinPhi1 = std::sin(phi1);
1024 cosPhi1 = std::cos(phi1); 1024 cosPhi1 = std::cos(phi1);
1025 G4double phi2 = CLHEP::twopi*G4UniformRand( 1025 G4double phi2 = CLHEP::twopi*G4UniformRand();
1026 sinPhi2 = std::sin(phi2); 1026 sinPhi2 = std::sin(phi2);
1027 cosPhi2 = std::cos(phi2); 1027 cosPhi2 = std::cos(phi2);
1028 1028
1029 // compute final direction realtive to z-di 1029 // compute final direction realtive to z-dir
1030 u2 = sinTheta2*cosPhi2; 1030 u2 = sinTheta2*cosPhi2;
1031 v2 = sinTheta2*sinPhi2; 1031 v2 = sinTheta2*sinPhi2;
1032 G4double u2p = cosTheta1*u2 + sinTheta1*cos 1032 G4double u2p = cosTheta1*u2 + sinTheta1*cosTheta2;
1033 uss = u2p*cosPhi1 - v2*sinPhi1; 1033 uss = u2p*cosPhi1 - v2*sinPhi1;
1034 vss = u2p*sinPhi1 + v2*cosPhi1; 1034 vss = u2p*sinPhi1 + v2*cosPhi1;
1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2; 1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2;
1036 1036
1037 // set new direction (is scattering frame) 1037 // set new direction (is scattering frame)
1038 fTheNewDirection.set(uss,vss,wss); 1038 fTheNewDirection.set(uss,vss,wss);
1039 1039
1040 // set the fTheZPathLenght if we don't samp 1040 // set the fTheZPathLenght if we don't sample displacement and
1041 // we should do everything at the step-limi 1041 // we should do everything at the step-limit-phase before we return
1042 if(fIsNoDisplace && fIsEverythingWasDone) 1042 if(fIsNoDisplace && fIsEverythingWasDone)
1043 fTheZPathLenght = fTheTrueStepLenght; 1043 fTheZPathLenght = fTheTrueStepLenght;
1044 1044
1045 // in optimized-mode if the current-safety 1045 // in optimized-mode if the current-safety > current-range we do not use dispalcement
1046 if(fIsNoDisplace) 1046 if(fIsNoDisplace)
1047 return; 1047 return;
1048 1048
1049 /////////////////////////////////////////// 1049 //////////////////////////////////////////////////////////////////////
1050 // Compute final position 1050 // Compute final position
1051 Qn1 *= fMCtoQ1; 1051 Qn1 *= fMCtoQ1;
1052 if (gIsUseAccurate) { 1052 if (gIsUseAccurate) {
1053 // correction parameter 1053 // correction parameter
1054 G4double par =1.; 1054 G4double par =1.;
1055 if(Qn1<0.7) par = 1.; 1055 if(Qn1<0.7) par = 1.;
1056 else if (Qn1<7.0) par = -0.031376*Qn1+1. 1056 else if (Qn1<7.0) par = -0.031376*Qn1+1.01356;
1057 else par = 0.79; 1057 else par = 0.79;
1058 1058
1059 // Moments with energy loss correction 1059 // Moments with energy loss correction
1060 // --first the uncorrected (for energy l 1060 // --first the uncorrected (for energy loss) values of gamma, eta, a1=a2=0.5*(1-eta), delta
1061 // gamma = G_2/G_1 based on G2 computed 1061 // gamma = G_2/G_1 based on G2 computed from A by using the Wentzel DCS form of G2
1062 G4double loga = G4Log(1.0+1.0/fScrA); 1062 G4double loga = G4Log(1.0+1.0/fScrA);
1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA 1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA)*(loga*(1.0 + 2.0*fScrA) - 2.0)/fG1;
1064 gamma *= fMCtoG2PerG1; 1064 gamma *= fMCtoG2PerG1;
1065 // sample eta from p(eta)=2*eta i.e. P(e 1065 // sample eta from p(eta)=2*eta i.e. P(eta) = eta_square ;-> P(eta) = rand --> eta = sqrt(rand)
1066 G4double eta = std::sqrt(G4UniformRan 1066 G4double eta = std::sqrt(G4UniformRand());
1067 G4double eta1 = 0.5*(1 - eta); // use 1067 G4double eta1 = 0.5*(1 - eta); // used more than once
1068 // 0.5 +sqrt(6)/6 = 0.9082483; 1068 // 0.5 +sqrt(6)/6 = 0.9082483;
1069 // 1/(4*sqrt(6)) = 0.1020621; 1069 // 1/(4*sqrt(6)) = 0.1020621;
1070 // (4-sqrt(6)/(24*sqrt(6))) = 0.02637471 1070 // (4-sqrt(6)/(24*sqrt(6))) = 0.026374715
1071 // delta = 0.9082483-(0.1020621-0.026374 1071 // delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1 without energy loss cor.
1072 G4double delta = 0.9082483-(0.1020621-0 1072 G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1;
1073 1073
1074 // compute alpha1 and alpha2 for energy 1074 // compute alpha1 and alpha2 for energy loss correction
1075 G4double temp1 = 2.0 + tau; 1075 G4double temp1 = 2.0 + tau;
1076 G4double temp = (2.0+tau*temp1)/((tau+1 1076 G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1);
1077 //Take logarithmic dependence 1077 //Take logarithmic dependence
1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga 1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0));
1079 temp = temp * epsm; 1079 temp = temp * epsm;
1080 temp1 = 1.0 - temp; 1080 temp1 = 1.0 - temp;
1081 delta = delta + 0.40824829*(eps0*(tau+1. 1081 delta = delta + 0.40824829*(eps0*(tau+1.0)/((tau+2.0)*
1082 (loga*(1.0+fScrA)-1.0)*(loga*(1. 1082 (loga*(1.0+fScrA)-1.0)*(loga*(1.0+2.0*fScrA)-2.0)) - 0.25*temp*temp);
1083 G4double b = eta*delta; 1083 G4double b = eta*delta;
1084 G4double c = eta*(1.0-delta); 1084 G4double c = eta*(1.0-delta);
1085 1085
1086 //Calculate transport direction cosines: 1086 //Calculate transport direction cosines:
1087 // ut,vt,wt is the final position divide 1087 // ut,vt,wt is the final position divided by the true step length
1088 G4double w1v2 = cosTheta1*v2; 1088 G4double w1v2 = cosTheta1*v2;
1089 G4double ut = b*sinTheta1*cosPhi1 + c* 1089 G4double ut = b*sinTheta1*cosPhi1 + c*(cosPhi1*u2 - sinPhi1*w1v2) + eta1*uss*temp1;
1090 G4double vt = b*sinTheta1*sinPhi1 + c* 1090 G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vss*temp1;
1091 G4double wt = eta1*(1+temp) + b* 1091 G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*wss*temp1;
1092 1092
1093 // long step correction 1093 // long step correction
1094 ut *=par; 1094 ut *=par;
1095 vt *=par; 1095 vt *=par;
1096 wt *=par; 1096 wt *=par;
1097 1097
1098 // final position relative to the pre-st 1098 // final position relative to the pre-step point in the scattering frame
1099 // ut = x_f/s so needs to multiply by s 1099 // ut = x_f/s so needs to multiply by s
1100 x_coord = ut*fTheTrueStepLenght; 1100 x_coord = ut*fTheTrueStepLenght;
1101 y_coord = vt*fTheTrueStepLenght; 1101 y_coord = vt*fTheTrueStepLenght;
1102 z_coord = wt*fTheTrueStepLenght; 1102 z_coord = wt*fTheTrueStepLenght;
1103 1103
1104 if(fIsEverythingWasDone){ 1104 if(fIsEverythingWasDone){
1105 // We sample in the step limit so set 1105 // We sample in the step limit so set fTheZPathLenght = transportDistance
1106 // and lateral displacement (x_coord,y 1106 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1107 //Calculate transport distance 1107 //Calculate transport distance
1108 G4double transportDistance = std::sqr 1108 G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
1109 // protection 1109 // protection
1110 if(transportDistance>fTheTrueStepLengh 1110 if(transportDistance>fTheTrueStepLenght)
1111 transportDistance = fTheTrueStepLen 1111 transportDistance = fTheTrueStepLenght;
1112 fTheZPathLenght = transportDistance; 1112 fTheZPathLenght = transportDistance;
1113 } 1113 }
1114 // else:: we sample in the DoIt so 1114 // else:: we sample in the DoIt so
1115 // the fTheZPathLenght was already 1115 // the fTheZPathLenght was already set and was taken as transport along zet
1116 fTheDisplacementVector.set(x_coord,y_coo 1116 fTheDisplacementVector.set(x_coord,y_coord,z_coord-fTheZPathLenght);
1117 } else { 1117 } else {
1118 // compute zz = <z>/tPathLength 1118 // compute zz = <z>/tPathLength
1119 // s -> true-path-length 1119 // s -> true-path-length
1120 // z -> geom-path-length:: when PRESTA i 1120 // z -> geom-path-length:: when PRESTA is used z =(def.) <z>
1121 // r -> lateral displacement = s/2 sin(t 1121 // r -> lateral displacement = s/2 sin(theta) => x_f = r cos(phi); y_f = r sin(phi)
1122 G4double zz = 0.0; 1122 G4double zz = 0.0;
1123 if(fIsEverythingWasDone){ 1123 if(fIsEverythingWasDone){
1124 // We sample in the step limit so set 1124 // We sample in the step limit so set fTheZPathLenght = transportDistance
1125 // and lateral displacement (x_coord, 1125 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1126 if(Qn1<0.1) { // use 3-order Taylor a 1126 if(Qn1<0.1) { // use 3-order Taylor approximation of (1-exp(-x))/x around x=0
1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666 1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1)); // 1/6 =0.166..7 ; 1/24=0.041..
1128 } else { 1128 } else {
1129 zz = (1.-G4Exp(-Qn1))/Qn1; 1129 zz = (1.-G4Exp(-Qn1))/Qn1;
1130 } 1130 }
1131 } else { 1131 } else {
1132 // we sample in the DoIt so 1132 // we sample in the DoIt so
1133 // the fTheZPathLenght was already se 1133 // the fTheZPathLenght was already set and was taken as transport along zet
1134 zz = fTheZPathLenght/fTheTrueStepLeng 1134 zz = fTheZPathLenght/fTheTrueStepLenght;
1135 } 1135 }
1136 1136
1137 G4double rr = (1.-zz*zz)/(1.-wss*wss); / 1137 G4double rr = (1.-zz*zz)/(1.-wss*wss); // s^2 >= <z>^2+r^2 :: where r^2 = s^2/4 sin^2(theta)
1138 if(rr >= 0.25) rr = 0.25; // 1138 if(rr >= 0.25) rr = 0.25; // (1-<z>^2/s^2)/sin^2(theta) >= r^2/(s^2 sin^2(theta)) = 1/4 must hold
1139 G4double rperp = fTheTrueStepLenght*std: 1139 G4double rperp = fTheTrueStepLenght*std::sqrt(rr); // this is r/sint
1140 x_coord = rperp*uss; 1140 x_coord = rperp*uss;
1141 y_coord = rperp*vss; 1141 y_coord = rperp*vss;
1142 z_coord = zz*fTheTrueStepLenght; 1142 z_coord = zz*fTheTrueStepLenght;
1143 1143
1144 if(fIsEverythingWasDone){ 1144 if(fIsEverythingWasDone){
1145 G4double transportDistance = std::sqrt 1145 G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord);
1146 fTheZPathLenght = transportDistance; 1146 fTheZPathLenght = transportDistance;
1147 } 1147 }
1148 1148
1149 fTheDisplacementVector.set(x_coord,y_coo 1149 fTheDisplacementVector.set(x_coord,y_coord,z_coord- fTheZPathLenght);
1150 } 1150 }
1151 } 1151 }
1152 1152