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