| 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 93663 2015-10-28 09:50:49Z 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 tlimit = 1.e+10*mm;
174 presafety = 0.*mm; << 136 presafety = 0.*mm;
175 // << 137
176 particle = nullptr; << 138 particle = 0;
177 theManager = G4LossTableManager: << 139 theManager = G4LossTableManager::Instance();
178 firstStep = true; << 140 firstStep = true;
179 currentKinEnergy = 0.0; << 141
180 currentRange = 0.0; << 142 tlimitminfix2 = 1.*nm;
181 // << 143 par1=par2=par3= 0.0;
182 tlimitminfix2 = 1.*nm; << 144 tausmall = 1.e-16;
183 tausmall = 1.e-16; << 145 mass = electron_mass_c2;
184 mass = electron_mass_c2; << 146 taulim = 1.e-6;
185 taulim = 1.e-6; << 147
186 // << 148 currentCouple = 0;
187 currentCouple = nullptr; << 149 fParticleChange = 0;
188 fParticleChange = nullptr; << 150
189 // << 151 // by default Moliere screeing parameter will be used but it can be set to the
190 fZeff = 1.; << 152 // PWA screening before calling the Initialise method
191 // << 153 fIsUsePWATotalXsecData = false;
192 par1 = 0.; << 154
193 par2 = 0.; << 155 //*****************
194 par3 = 0.; << 156 fLambda0 = 0.0; // elastic mean free path
195 // << 157 fLambda1 = 0.0; // first transport mean free path
196 // Moliere screeing parameter will be used a << 158 fScrA = 0.0; // screening parameter
197 // appalied to the integrated quantities (sc << 159 fG1 = 0.0; // first transport coef.
198 // and second moments) derived from the corr << 160 //******************
199 // this PWA correction is ignored if Mott-co << 161 fTheTrueStepLenght = 0.;
200 // Mott-correction contains all these correc << 162 fTheTransportDistance = 0.;
201 fIsUsePWACorrection = true; << 163 fTheZPathLenght = 0.;
202 // <<
203 fIsUseMottCorrection = false; <<
204 // <<
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.); 164 fTheDisplacementVector.set(0.,0.,0.);
219 fTheNewDirection.set(0.,0.,1.); 165 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 166
>> 167 fIsEverythingWasDone = false;
>> 168 fIsMultipleSacettring = false;
>> 169 fIsSingleScattering = false;
>> 170 fIsEndedUpOnBoundary = false;
>> 171 fIsNoScatteringInMSC = false;
>> 172 fIsNoDisplace = false;
>> 173
>> 174 fZeff = 1.;
>> 175 //+++++++++++++
>> 176
>> 177 rndmEngineMod = G4Random::getTheEngine();
>> 178 }
236 179
237 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaund << 180 ////////////////////////////////////////////////////////////////////////////////
238 if (IsMaster()) { << 181 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel(){
239 if (fGSTable) { << 182 if(IsMaster()){
240 delete fGSTable; << 183 if(fgGSTable){
241 fGSTable = nullptr; << 184 delete fgGSTable;
>> 185 fgGSTable = 0;
242 } 186 }
243 if (fPWACorrection) { << 187 if(fgPWAXsecTable){
244 delete fPWACorrection; << 188 delete fgPWAXsecTable;
245 fPWACorrection = nullptr; << 189 fgPWAXsecTable = 0;
246 } 190 }
247 } 191 }
248 } 192 }
249 193
250 194
251 void G4GoudsmitSaundersonMscModel::Initialise( << 195 ////////////////////////////////////////////////////////////////////////////////
>> 196 void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
>> 197 const G4DataVector&){
>> 198 //skindepth=skin*stepmin;
252 SetParticle(p); 199 SetParticle(p);
253 InitialiseParameters(p); <<
254 // -create GoudsmitSaundersonTable and init <<
255 // Mott-correction was required <<
256 if (IsMaster()) { <<
257 // get the Mott-correction flag from EmPar <<
258 if (G4EmParameters::Instance()->UseMottCor <<
259 fIsUseMottCorrection = true; <<
260 } <<
261 // Mott-correction includes other way of P <<
262 // when Mott-correction is activated by th <<
263 if (fIsUseMottCorrection) { <<
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 } <<
276 // create GS-table <<
277 G4bool isElectron = true; <<
278 if (p->GetPDGCharge()>0.) { <<
279 isElectron = false; <<
280 } <<
281 fGSTable = new G4GoudsmitSaundersonTable(i <<
282 // G4GSTable will be initialised: <<
283 // - Screened-Rutherford DCS based GS angu <<
284 // - Mott-correction will be initialised i <<
285 fGSTable->SetOptionMottCorrection(fIsUseMo <<
286 // - set PWA correction (correction to int <<
287 fGSTable->SetOptionPWACorrection(fIsUsePWA <<
288 // init <<
289 fGSTable->Initialise(LowEnergyLimit(),High <<
290 // create PWA corrections table if it was <<
291 if (fIsUsePWACorrection) { <<
292 fPWACorrection = new G4GSPWACorrections( <<
293 fPWACorrection->Initialise(); <<
294 } <<
295 } <<
296 fParticleChange = GetParticleChangeForMSC(p) 200 fParticleChange = GetParticleChangeForMSC(p);
297 } <<
298 201
>> 202 //
>> 203 // -create static GoudsmitSaundersonTable and PWATotalXsecTable is necessary
>> 204 //
>> 205 if(IsMaster()) {
299 206
300 void G4GoudsmitSaundersonMscModel::InitialiseL << 207 // Could delete as well if they are exist but then the same data are reloaded
301 fGSTable = static_cast<G4Goud << 208 // in case of e+ for example.
302 fIsUseMottCorrection = static_cast<G4Goud << 209 /* if(fgGSTable) {
303 fIsUsePWACorrection = static_cast<G4Goud << 210 delete fgGSTable;
304 fPWACorrection = static_cast<G4Goud << 211 fgGSTable = 0;
305 } << 212 }
>> 213 if(fgPWAXsecTable) {
>> 214 delete fgPWAXsecTable;
>> 215 fgPWAXsecTable = 0;
>> 216 }
>> 217 */
>> 218 if(!fgGSTable)
>> 219 fgGSTable = new G4GoudsmitSaundersonTable();
306 220
>> 221 // G4GSTable will be initialised (data are loaded) only if they are not there yet.
>> 222 fgGSTable->Initialise();
307 223
308 // computes macroscopic first transport cross << 224 if(fIsUsePWATotalXsecData) {
309 G4double G4GoudsmitSaundersonMscModel::CrossSe << 225 if(!fgPWAXsecTable)
310 const << 226 fgPWAXsecTable = new G4PWATotalXsecTable();
311 G4dou << 227
312 G4dou << 228 // Makes sure that all (and only those) atomic PWA data are in memory that
313 G4dou << 229 // are in the current MaterilaTable.
314 G4double xsecTr1 = 0.; // cross section per << 230 fgPWAXsecTable->Initialise();
315 // << 231 }
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 } 232 }
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 } 233 }
361 234
362 235
>> 236 ///////////////////////////////////////////////////////////////////////////////
363 // gives back the first transport mean free pa 237 // gives back the first transport mean free path in internal G4 units
364 G4double 238 G4double
365 G4GoudsmitSaundersonMscModel::GetTransportMean 239 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePath(const G4ParticleDefinition* /*partdef*/,
366 240 G4double kineticEnergy) {
367 // kinetic energy is assumed to be in Geant4 241 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
368 G4double efEnergy = kineticEnergy; 242 G4double efEnergy = kineticEnergy;
369 // << 243 if(efEnergy<lowKEnergy) efEnergy=lowKEnergy;
>> 244 if(efEnergy>highKEnergy)efEnergy=highKEnergy;
>> 245
>> 246 ///////////////////////////////////////////
370 const G4Material* mat = currentCouple->GetM 247 const G4Material* mat = currentCouple->GetMaterial();
371 // << 248
372 fLambda0 = 0.0; // elastic mean free path 249 fLambda0 = 0.0; // elastic mean free path
373 fLambda1 = 0.0; // first transport mean free 250 fLambda1 = 0.0; // first transport mean free path
374 fScrA = 0.0; // screening parameter 251 fScrA = 0.0; // screening parameter
375 fG1 = 0.0; // first transport coef. 252 fG1 = 0.0; // first transport coef.
376 253
377 // use Moliere's screening (with Mott-corret << 254 if(fIsUsePWATotalXsecData) {
378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* << 255 // use PWA total xsec data to determine screening and lambda0 lambda1
379 // total mometum square << 256 const G4ElementVector* theElementVector = mat->GetElementVector();
380 G4double pt2 = efEnergy*(efEnergy+2.0*el << 257 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
381 // beta square << 258 G4int nelm = mat->GetNumberOfElements();
382 G4double beta2 = pt2/(pt2+electron_mass_c2 << 259
383 // current material index << 260 int elowindx = fgPWAXsecTable->GetPWATotalXsecForZet( G4lrint((*theElementVector)[0]->GetZ()) )
384 G4int matindx = (G4int)mat->GetIndex(); << 261 ->GetPWATotalXsecEnergyBinIndex(efEnergy);
385 // Moliere's b_c << 262 for(G4int i=0;i<nelm;i++){
386 G4double bc = fGSTable->GetMoliereBc(ma << 263 int zet = G4lrint((*theElementVector)[i]->GetZ());
387 // get the Mott-correcton factors if Mott-co << 264 // total elastic scattering cross section in Geant4 internal length2 units
388 fMCtoScrA = 1.0; << 265 int indx = (G4int)(1.5 + charge*1.5); // specify that want to get the total elastic scattering xsection
389 fMCtoQ1 = 1.0; << 266 double elxsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
390 fMCtoG2PerG1 = 1.0; << 267 // first transport cross section in Geant4 internal length2 units
391 G4double scpCor = 1.0; << 268 indx = (G4int)(2.5 + charge*1.5); // specify that want to get the first tarnsport xsection
392 if (fIsUseMottCorrection) { << 269 double t1xsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
393 fGSTable->GetMottCorrectionFactors(G4Log(e << 270 fLambda0 += (theAtomNumDensityVector[i]*elxsec);
394 scpCor = fGSTable->ComputeScatteringPowerC << 271 fLambda1 += (theAtomNumDensityVector[i]*t1xsec);
395 } else if (fIsUsePWACorrection) { << 272 }
396 fPWACorrection->GetPWACorrectionFactors(G4 << 273 if(fLambda0>0.0) { fLambda0 =1./fLambda0; }
397 // scpCor = fGSTable->ComputeScatteringPow << 274 if(fLambda1>0.0) { fLambda1 =1./fLambda1; }
>> 275
>> 276 // determine screening parameter such that it gives back PWA G1
>> 277 fG1=0.0;
>> 278 if(fLambda1>0.0) { fG1 = fLambda0/fLambda1; }
>> 279
>> 280 fScrA = fgGSTable->GetScreeningParam(fG1);
>> 281 } else {
>> 282 // Get SCREENING FROM MOLIER
>> 283 // total mometum square in Geant4 internal energy2 units which is MeV2
>> 284
>> 285 // below 1 keV it can give bananas so prevet (1 keV)
>> 286 if(efEnergy<0.001) efEnergy = 0.001;
>> 287
>> 288 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
>> 289 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
>> 290 G4int matindx = mat->GetIndex();
>> 291 G4double bc = fgGSTable->GetMoliereBc(matindx);
>> 292 fScrA = fgGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc);
>> 293 // total elastic mean free path in Geant4 internal lenght units
>> 294 fLambda0 = beta2/bc/(1.+fScrA);
>> 295 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/fScrA+1.0)-1.0);
>> 296 fLambda1 = fLambda0/fG1;
398 } 297 }
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 298
415 return fLambda1; 299 return fLambda1;
416 } 300 }
417 301
418 << 302 ////////////////////////////////////////////////////////////////////////////////
419 G4double 303 G4double
420 G4GoudsmitSaundersonMscModel::GetTransportMean 304 G4GoudsmitSaundersonMscModel::GetTransportMeanFreePathOnly(const G4ParticleDefinition* /*partdef*/,
421 305 G4double kineticEnergy) {
422 // kinetic energy is assumed to be in Geant4 306 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
423 G4double efEnergy = kineticEnergy; 307 G4double efEnergy = kineticEnergy;
424 // << 308 if(efEnergy<lowKEnergy) efEnergy=lowKEnergy;
>> 309 if(efEnergy>highKEnergy)efEnergy=highKEnergy;
>> 310
>> 311 ///////////////////////////////////////////
425 const G4Material* mat = currentCouple->GetM 312 const G4Material* mat = currentCouple->GetMaterial();
426 // << 313
427 G4double lambda0 = 0.0; // elastc mean free 314 G4double lambda0 = 0.0; // elastc mean free path
428 G4double lambda1 = 0.0; // first transport m 315 G4double lambda1 = 0.0; // first transport mean free path
429 G4double scrA = 0.0; // screening paramet 316 G4double scrA = 0.0; // screening parametr
430 G4double g1 = 0.0; // first transport m 317 G4double g1 = 0.0; // first transport mean free path
431 318
432 // use Moliere's screening (with Mott-corret << 319 if(fIsUsePWATotalXsecData) {
433 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* << 320 // use PWA total xsec data to determine screening and lambda0 lambda1
434 // total mometum square in Geant4 internal e << 321 const G4ElementVector* theElementVector = mat->GetElementVector();
435 G4double pt2 = efEnergy*(efEnergy+2.0*el << 322 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
436 G4double beta2 = pt2/(pt2+electron_mass_c2 << 323 G4int nelm = mat->GetNumberOfElements();
437 G4int matindx = (G4int)mat->GetIndex(); << 324
438 G4double bc = fGSTable->GetMoliereBc(ma << 325 int elowindx = fgPWAXsecTable->GetPWATotalXsecForZet( G4lrint((*theElementVector)[0]->GetZ()) )
439 // get the Mott-correcton factors if Mott-co << 326 ->GetPWATotalXsecEnergyBinIndex(efEnergy);
440 G4double mctoScrA = 1.0; << 327 for(G4int i=0;i<nelm;i++){
441 G4double mctoQ1 = 1.0; << 328 int zet = G4lrint((*theElementVector)[i]->GetZ());
442 G4double mctoG2PerG1 = 1.0; << 329 // total elastic scattering cross section in Geant4 internal length2 units
443 G4double scpCor = 1.0; << 330 // int indx = (G4int)(1.5 + charge*1.5); // specify that want to get the total elastic scattering xsection
444 if (fIsUseMottCorrection) { << 331 // double elxsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
445 fGSTable->GetMottCorrectionFactors(G4Log(e << 332 // first transport cross section in Geant4 internal length2 units
446 scpCor = fGSTable->ComputeScatteringPowerC << 333 int indx = (G4int)(2.5 + charge*1.5); // specify that want to get the first tarnsport xsection
447 } else if (fIsUsePWACorrection) { << 334 double t1xsec = (fgPWAXsecTable->GetPWATotalXsecForZet(zet))->GetInterpXsec(efEnergy,elowindx,indx);
448 fPWACorrection->GetPWACorrectionFactors(G4 << 335 // fLambda0 += (theAtomNumDensityVector[i]*elxsec);
449 // scpCor = fGSTable->ComputeScatteringPow << 336 lambda1 += (theAtomNumDensityVector[i]*t1xsec);
>> 337 }
>> 338 // if(fLambda0>0.0) { fLambda0 =1./fLambda0; }
>> 339 if(lambda1>0.0) { lambda1 =1./lambda1; }
>> 340
>> 341 } else {
>> 342 // Get SCREENING FROM MOLIER
>> 343 // total mometum square in Geant4 internal energy2 units which is MeV2
>> 344
>> 345 // below 1 keV it can give bananas so prevet (1 keV)
>> 346 if(efEnergy<0.001) efEnergy = 0.001;
>> 347
>> 348 G4double pt2 = efEnergy*(efEnergy+2.0*electron_mass_c2);
>> 349 G4double beta2 = pt2/(pt2+electron_mass_c2*electron_mass_c2);
>> 350 G4int matindx = mat->GetIndex();
>> 351 G4double bc = fgGSTable->GetMoliereBc(matindx);
>> 352 scrA = fgGSTable->GetMoliereXc2(matindx)/(4.0*pt2*bc);
>> 353 // total elastic mean free path in Geant4 internal lenght units
>> 354 lambda0 = beta2/bc/(1.+scrA);
>> 355 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scrA+1.0)-1.0);
>> 356 lambda1 = lambda0/g1;
450 } 357 }
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 358
457 return lambda1; 359 return lambda1;
458 } 360 }
459 361
460 << 362 ////////////////////////////////////////////////////////////////////////////////
461 void G4GoudsmitSaundersonMscModel::StartTracki << 363 void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track)
>> 364 {
462 SetParticle(track->GetDynamicParticle()->Get 365 SetParticle(track->GetDynamicParticle()->GetDefinition());
463 firstStep = true; 366 firstStep = true;
464 tlimit = tgeom = rangeinit = geombig; <<
465 rangeinit = 1.e+21; <<
466 } 367 }
467 368
468 369
469 G4double G4GoudsmitSaundersonMscModel::Compute << 370 ////////////////////////////////////////////////////////////////////////////////
470 << 371 G4double G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(
>> 372 const G4Track& track,
>> 373 G4double& currentMinimalStep){
>> 374 // some constant
>> 375 const G4double almostZero = 1.e-8;
>> 376 const G4double almostOne = 1.-almostZero;
>> 377 const G4double mlogAlmostZero = -1.*G4Log(almostZero);
>> 378 const G4double onemExpmOne = 1.-G4Exp(-1.0);
>> 379
471 G4double skindepth = 0.; 380 G4double skindepth = 0.;
472 // << 381
473 const G4DynamicParticle* dp = track.GetDynam 382 const G4DynamicParticle* dp = track.GetDynamicParticle();
474 G4StepPoint* sp = track.GetStep( << 383
475 G4StepStatus stepStatus = sp->GetStepSta << 384 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
476 currentCouple = track.GetMater << 385 G4StepStatus stepStatus = sp->GetStepStatus();
>> 386 currentCouple = track.GetMaterialCutsCouple();
477 SetCurrentCouple(currentCouple); 387 SetCurrentCouple(currentCouple);
478 currentMaterialIndex = (G4int)current << 388 currentMaterialIndex = currentCouple->GetIndex();
479 currentKinEnergy = dp->GetKinetic << 389 currentKinEnergy = dp->GetKineticEnergy();
480 currentRange = GetRange(parti << 390 currentRange = GetRange(particle,currentKinEnergy,currentCouple);
481 dp->G << 391
482 // elastic and first transport mfp, screenin << 392
483 // (Mott-correction will be used if it was r << 393
>> 394 // *** Elastic and first transport mfp, fScrA and fG1 are also set in this !!!
484 fLambda1 = GetTransportMeanFreePath(particle 395 fLambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
>> 396
>> 397 //
485 // Set initial values: 398 // Set initial values:
486 // : lengths are initialised to currentMini 399 // : lengths are initialised to currentMinimalStep which is the true, minimum
487 // step length from all other physics 400 // step length from all other physics
488 fTheTrueStepLenght = currentMinimalStep; 401 fTheTrueStepLenght = currentMinimalStep;
489 fTheTransportDistance = currentMinimalStep; 402 fTheTransportDistance = currentMinimalStep;
490 fTheZPathLenght = currentMinimalStep; 403 fTheZPathLenght = currentMinimalStep; // will need to be converted
491 fTheDisplacementVector.set(0.,0.,0.); 404 fTheDisplacementVector.set(0.,0.,0.);
492 fTheNewDirection.set(0.,0.,1.); 405 fTheNewDirection.set(0.,0.,1.);
493 406
494 // Can everything be done in the step limit 407 // Can everything be done in the step limit phase ?
495 fIsEverythingWasDone = false; 408 fIsEverythingWasDone = false;
496 // Multiple scattering needs to be sample ? 409 // Multiple scattering needs to be sample ?
497 fIsMultipleSacettring = false; 410 fIsMultipleSacettring = false;
498 // Single scattering needs to be sample ? 411 // Single scattering needs to be sample ?
499 fIsSingleScattering = false; 412 fIsSingleScattering = false;
500 // Was zero deflection in multiple scatterin 413 // Was zero deflection in multiple scattering sampling ?
501 fIsNoScatteringInMSC = false; 414 fIsNoScatteringInMSC = false;
502 // Do not care about displacement in MSC sam 415 // Do not care about displacement in MSC sampling
503 // ( used only in the case of gIsOptimizatio << 416 // ( used only in the case of fgIsOptimizationOn = true)
504 fIsNoDisplace = false; 417 fIsNoDisplace = false;
>> 418
505 // get pre-step point safety 419 // get pre-step point safety
506 presafety = sp->GetSafety(); 420 presafety = sp->GetSafety();
507 // << 421
508 fZeff = currentCouple->GetMaterial()->GetIon << 422 fZeff = currentCouple->GetMaterial()->GetTotNbOfElectPerVolume()/
509 // distance will take into account max-fluct << 423 currentCouple->GetMaterial()->GetTotNbOfAtomsPerVolume();
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));
512 // << 426
513 // Possible optimization : if the distance i << 427 // In G4VMscModel there is a member G4MscStepLimitType steppingAlgoritm;
514 // particle will never leave this volume -> << 428 if( steppingAlgorithm == fUseDistanceToBoundary || steppingAlgorithm == fUseSafety ) {
515 // as the effect of multiple elastic scatter << 429 // Possible optimization : if the range is samller than the safety -> the
516 // Important : this optimization can cause p << 430 // particle will never leave this volume -> dispalcement and angular
517 // in a bigger volume since MSC won't be don << 431 // deflection as the effects of multiple elastic scattering can be skipped
518 // distance < safety so don't use optimized- << 432 // Only true->geometrical->true pathlength conversion is done based on mean
519 if (gIsOptimizationOn && (distance<presafety << 433 // values. In order to be able to do this, only the maximum geometrical step
520 // Indicate that we need to do MSC after << 434 // length will be checked.
521 fIsMultipleSacettring = true; << 435 // Important : this optimization can cause problems if one does scoring
522 fIsNoDisplace = true; << 436 // in a bigger volume since MSC won't be done deep inside the volume when
523 } else if (steppingAlgorithm==fUseDistanceTo << 437 // range < safety so don't use optimized-mode in such case.
524 //Compute geomlimit (and presafety) : << 438 if( fgIsOptimizationOn && (distance < presafety) ) {
525 // - geomlimit will be: << 439 // Do path length conversion using the means and do NOT sample dispalcement.
526 // == the straight line distance to the << 440 // :: according to optimized-mode:
527 // longer than that << 441 // restrict the true step length such that the corresponding
528 // == a big value [geombig = 1.e50*mm] << 442 // transport length is not longer than the first transport mean free path:
529 // the straight line distance to the << 443 // -this is important for the g->t conversion (i.e. fTheZPathLenght ->
530 // - presafety will be updated as well << 444 // fTheTrueStepLenght) that will be invoked later based on the mean
531 // So the particle can travell 'gemlimit' << 445 // values. However, it leads to fTheTrueStepLenght >> lambda1 and
532 // line!) in its current direction: << 446 // the MSC model as condensed simulation model is valid about
533 // (1) before reaching a boundary (geomli << 447 // fTheTrueStepLenght ~ lambda1 that corresponds to <cos(theta)> ~ 1 rad
534 // (2) before reaching its current range << 448 // in one sub-step if we divide the true step length into two parts when
535 geomlimit = ComputeGeomLimit(track, presaf << 449 // sampling the angular distribution.
536 // Record that the particle is on a bounda << 450 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*mlogAlmostZero);
537 if ( (stepStatus==fGeomBoundary) || (stepS << 451 fTheZPathLenght = fLambda1*(1.-G4Exp(-fTheTrueStepLenght/fLambda1));
538 fIsWasOnBoundary = true; << 452 // Indicate that we need to do MSC after transportation and no dispalcement.
539 } << 453 fIsMultipleSacettring = true;
540 // Set skin depth = skin x elastic_mean_fr << 454 fIsNoDisplace = true;
541 skindepth = skin*fLambda0; << 455 } else {
542 // Init the flag that indicates that the p << 456
543 // distance from a boundary << 457
544 fIsInsideSkin = false; << 458 if(steppingAlgorithm == fUseDistanceToBoundary) {
545 // Check if we can try Single Scattering b << 459 geomlimit = ComputeGeomLimit(track, presafety, currentRange);
546 // distance from/to a boundary OR the curr <<
547 // shorter than skindepth. NOTICE: the lat <<
548 // because the MSC angular sampling is fin <<
549 // faster to try single scattering in case <<
550 if ((stepStatus==fGeomBoundary) || (presaf <<
551 // check if we are within skindepth dist <<
552 if ((stepStatus == fGeomBoundary) || (pr <<
553 fIsInsideSkin = true; <<
554 fIsWasOnBoundary = true; <<
555 } <<
556 //Try single scattering: <<
557 // - sample distance to next single scat <<
558 // - compare to current minimum length <<
559 // == if sslimit is the shorter: <<
560 // - set the step length to ssl <<
561 // - indicate that single scatt <<
562 // == else : nothing to do <<
563 //- in both cases, the step length was v <<
564 // true path length are the same <<
565 G4double sslimit = -1.*fLambda0*G4Log(G4 <<
566 // compare to current minimum step lengt <<
567 if (sslimit<fTheTrueStepLenght) { <<
568 fTheTrueStepLenght = sslimit; <<
569 fIsSingleScattering = true; <<
570 } <<
571 // short step -> true step length equal <<
572 fTheZPathLenght = fTheTrueStepLeng <<
573 // Set taht everything is done in step-l <<
574 // We will check if we need to perform t <<
575 // sampling i.e. if single elastic scatt <<
576 fIsEverythingWasDone = true; <<
577 } else { <<
578 // After checking we know that we cannot <<
579 // need to make an MSC step <<
580 // Indicate that we need to make and MSC <<
581 // do it now i.e. if presafety>final_tru <<
582 // fIsEverythingWasDone = false which in <<
583 // MSC after transportation. <<
584 fIsMultipleSacettring = true; <<
585 // Init the first-real-step falg: it wil <<
586 // non-single scattering step in this vo <<
587 fIsFirstRealStep = false; <<
588 // If previously the partcile was on bou <<
589 // well. When it is not within skin anym <<
590 // so we make the first real MSC step wi <<
591 if (fIsWasOnBoundary && !fIsInsideSkin) <<
592 // reset the 'was on boundary' indicat <<
593 fIsWasOnBoundary = false; <<
594 fIsFirstRealStep = true; <<
595 } <<
596 // If this is the first-real msc step (t <<
597 // skin) or this is the first step with <<
598 // primary): <<
599 // - set the initial range that will b <<
600 // (only in this volume, because aft <<
601 // first-real MSC step we will reset <<
602 // - don't let the partcile to cross th <<
603 if (firstStep || fIsFirstRealStep || ran <<
604 rangeinit = currentRange; <<
605 // If geomlimit < geombig than the par <<
606 // along its initial direction before <<
607 // Otherwise we can be sure that the p <<
608 // before reaching the boundary along <<
609 // geometrical limit appalied. [Howeve <<
610 // first or the first-real MSC step. A <<
611 // MSC step the direction will change <<
612 // But we will try to end up within sk <<
613 // the actual value of geomlimit(See l <<
614 // boundary).] <<
615 if (geomlimit<geombig) { <<
616 // transfrom straight line distance <<
617 // length based on the mean values ( <<
618 // first-transport mean free path i. <<
619 if ((1.-geomlimit/fLambda1)> 0.) { <<
620 geomlimit = -fLambda1*G4Log(1.-geo <<
621 } <<
622 // the 2-different case that could l <<
623 if (firstStep) { <<
624 tgeom = 2.*geomlimit/facgeom; <<
625 } else { <<
626 tgeom = geomlimit/facgeom; <<
627 } <<
628 } else { 460 } else {
629 tgeom = geombig; << 461 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
>> 462 geomlimit = presafety;
>> 463 }
>> 464
>> 465 // Set limit from range factor even in non-optimized mode.
>> 466 // Important: At lower energies the energy loss process continuous limit
>> 467 // will converge to range.
>> 468 // In this case, we will limit the step to 'facrange x range' as long
>> 469 // as the range is higher than safety i.e. the particle can leave the
>> 470 // current volume. One can completely eliminate this limit by setting
>> 471 // 'facrange = 1.' .
>> 472 if(currentRange > presafety)
>> 473 fTheTrueStepLenght = std::min(facrange*currentRange, fTheTrueStepLenght);
>> 474
>> 475 // Compute skin-depth; fLambda0(elastic-mfp) was already set in
>> 476 // GetTransportMeanFreePath(...)
>> 477 skindepth = skin*fLambda0;
>> 478
>> 479 // Limit the true step -length to the first transport mean free path that
>> 480 // is the limit of the MSC model. Note: if the true step lenght is about
>> 481 // half transport mean free path i.e. t = 0.5 x lambda_1 the corresponding
>> 482 // mean angular deflection <cos(\theta)> = exp(-t/lambda_1) -> \theta ~
>> 483 // 0.919 radian (52.6 degree) that is the aplicability limit of any
>> 484 // condensed history techniques. The result will be more and more
>> 485 // questionable when condensed history techniques are applied for even
>> 486 // longer true step lengts.
>> 487 // if( (currentKinEnergy > 0.3) || !fgIsOptimizationOn)
>> 488 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*almostOne);
>> 489
>> 490 // Check if we are within skin depth to/from boundary or making a step
>> 491 // that is shorter than skindepth distance => try single scattering::
>> 492 // - if stepStatus is fGeomBoundary -> pre-step point is on boundary
>> 493 // - in case of initStep pre-step point status is fUndefined:
>> 494 // at the moment I always call the navigator CheckNextStep even from
>> 495 // the WORLD that will set presafety; if presafety = 0 then we are on
>> 496 // the boundary at initStep
>> 497 // - if the currentMinimalStep is already shorter than skindepth we do
>> 498 // single scattering. NOTICE: it has only efficieny reasons because
>> 499 // the angular sampling is fine for any short steps but much faster to
>> 500 // try single scattering in case of short steps.
>> 501 if( (stepStatus == fGeomBoundary) || (presafety < skindepth) || (fTheTrueStepLenght < skindepth)){
>> 502 //Try single scattering:
>> 503 // - sample distance to next single scattering interaction
>> 504 // - compare to current minimum length
>> 505 // * if ss-length is the shorter:
>> 506 // - set the step length to the ss-length
>> 507 // - indicate that single scattering needs to be done
>> 508 // * else : nothing to do
>> 509 G4double sslimit = -1.*fLambda0*G4Log(G4UniformRand());
>> 510 // compare to current minimum step length
>> 511 if(sslimit < fTheTrueStepLenght) {
>> 512 fTheTrueStepLenght = sslimit;
>> 513 fIsSingleScattering = true;
>> 514 // zPathLength = tPathLength
>> 515 // single scattering needs to be done
>> 516 } //else {
>> 517 //- tPathLength = tPathLength
>> 518 //- zPathLength = tPathLength
>> 519 //-> nothing to do
>> 520 //}
>> 521
>> 522 // short step -> true step length equal to geometrical path length
>> 523 fTheZPathLenght = fTheTrueStepLenght;
>> 524 //set taht everything is done in step-limit phase
>> 525 fIsEverythingWasDone = true;
>> 526 } else {
>> 527 // Compare the true step length to presafty to see if we can do safe MSC step
>> 528 // - if yes, then performe MSC here
>> 529 // - if no, postpone MSC sampling after transportation
>> 530 if(fTheTrueStepLenght < presafety) {
>> 531 // => we can do safe MSC step because for sure we don't reach boundary
>> 532 // do multiple scattering: MSC will give the corresponding transport
>> 533 // length i.e. fTheZPathLenght (and new direction, displacement)
>> 534 fIsMultipleSacettring = true;
>> 535 //set taht everything was done in step-limit phase
>> 536 fIsEverythingWasDone = true;
>> 537 } else {
>> 538 // => we cannot be sure at this point what will be the true step length
>> 539 // because transportation can hit boundary. The two step limit version
>> 540 // will be different but both case MSC will need to be done:
>> 541 // before(fUseSafety) OR after(fUseDistanceToBoundary) transportation.
>> 542 fIsMultipleSacettring = true;
>> 543 // Important that fIsEverythingWasDone remains false here;
>> 544
>> 545 // two versions will be different here
>> 546 if(geomlimit > presafety) {
>> 547 // => It can happen only if fUseDistanceToBoundary stepping is used
>> 548 // - limit geomlimit to a geometrical value, that corresponds
>> 549 // to the maximum true step length i.e. to lambda1
>> 550 // - convert the limited geomlimit to true length
>> 551 // - take minimum of this and the current minimal true step length
>> 552
>> 553 // if geom limit is not huge(real distance to boundary) then shorten it by almost skindepth
>> 554 // => we try to end up within the skindepth to the boundary and not on the boundary.
>> 555 if( geomlimit < geombig )
>> 556 geomlimit -= 0.999*skindepth;
>> 557
>> 558 // Set a maximum geometrical limit: a straight line distance that
>> 559 // corresponds to the maximal true step length i.e. lambda1 of
>> 560 // the MSC model (it is set based on the mean value).
>> 561 geomlimit = std::min(geomlimit, fLambda1*onemExpmOne);
>> 562 // convert to true-one
>> 563 geomlimit = -fLambda1*G4Log(1.-geomlimit/fLambda1);
>> 564 fTheTrueStepLenght = std::min(fTheTrueStepLenght,geomlimit);
>> 565 } else {
>> 566 // => It can happen only if geomlimit is the safety i.e. fUseSafety
>> 567 // steppingAlgorithm.
>> 568 // - if we are here, then presafety is for sure shorther than lambda1 so
>> 569 // we don't need to restrict
>> 570 // -limit the true step-length to presafety: for sure we can
>> 571 // perform MSC before transportation.
>> 572 fTheTrueStepLenght = presafety; // geomlimit = presafty
>> 573 fIsEverythingWasDone = true;
>> 574 }
>> 575 }
630 } 576 }
631 } 577 }
632 // True step length limit from range fac <<
633 // range is used that was set at the fir <<
634 // in this volume with this particle. <<
635 tlimit = facrange*rangeinit; <<
636 // Take the minimum of the true step len <<
637 // geometrical constraint or range-facto <<
638 tlimit = std::min(tlimit,tgeom); <<
639 // Step reduction close to boundary: we <<
640 // from the boundary ( Notice: in case o <<
641 // because geomlimit is the straigth lin <<
642 // the currect direction (if geomlimit<g <<
643 // change the initial direction. So te p <<
644 // before in a different direction. Howe <<
645 // path length to this (straight line) l <<
646 // transport distance (straight line) wi <<
647 // geomlimit-0.999*skindepth after the c <<
648 if (geomlimit<geombig) { <<
649 tlimit = std::min(tlimit, geomlimit-0. <<
650 } <<
651 // randomize 1st step or 1st 'normal' st <<
652 if (firstStep || fIsFirstRealStep) { <<
653 fTheTrueStepLenght = std::min(fTheTrue <<
654 } else { <<
655 fTheTrueStepLenght = std::min(fTheTrue <<
656 } <<
657 } <<
658 } else if (steppingAlgorithm==fUseSafetyPlus <<
659 presafety = ComputeSafety(sp->GetPosition <<
660 geomlimit = presafety; <<
661 // Set skin depth = skin x elastic_mean_fr <<
662 skindepth = skin*fLambda0; <<
663 // Check if we can try Single Scattering b <<
664 // distance from/to a boundary OR the curr <<
665 // shorter than skindepth. NOTICE: the lat <<
666 // because the MSC angular sampling is fin <<
667 // faster to try single scattering in case <<
668 if ((stepStatus==fGeomBoundary) || (presaf <<
669 //Try single scattering: <<
670 // - sample distance to next single scat <<
671 // - compare to current minimum length <<
672 // == if sslimit is the shorter: <<
673 // - set the step length to ssl <<
674 // - indicate that single scatt <<
675 // == else : nothing to do <<
676 //- in both cases, the step length was v <<
677 // true path length are the same <<
678 G4double sslimit = -1.*fLambda0*G4Log(G4 <<
679 // compare to current minimum step lengt <<
680 if (sslimit<fTheTrueStepLenght) { <<
681 fTheTrueStepLenght = sslimit; <<
682 fIsSingleScattering = true; <<
683 } <<
684 // short step -> true step length equal <<
685 fTheZPathLenght = fTheTrueStepLeng <<
686 // Set taht everything is done in step-l <<
687 // We will check if we need to perform t <<
688 // sampling i.e. if single elastic scatt <<
689 fIsEverythingWasDone = true; <<
690 } else { <<
691 // After checking we know that we cannot <<
692 // need to make an MSC step <<
693 // Indicate that we need to make and MSC <<
694 fIsMultipleSacettring = true; <<
695 fIsEverythingWasDone = true; <<
696 // limit from range factor <<
697 fTheTrueStepLenght = std::min(fTheTru <<
698 // never let the particle go further tha <<
699 // if we are here we are out of the skin <<
700 if (fTheTrueStepLenght>presafety) { <<
701 fTheTrueStepLenght = std::min(fTheTrue <<
702 } <<
703 // make sure that we are still within th <<
704 // i.e. true step length is not longer t <<
705 // We schould take into account energy l <<
706 // step length as well. So let it 0.5 x <<
707 fTheTrueStepLenght = std::min(fTheTrueSt <<
708 } <<
709 } else { 578 } else {
>> 579
710 // This is the default stepping algorithm: 580 // This is the default stepping algorithm: the fastest but the least
711 // accurate that corresponds to fUseSafety 581 // accurate that corresponds to fUseSafety in Urban model. Note, that GS
712 // model can handle any short steps so we 582 // model can handle any short steps so we do not need the minimum limits
713 // << 583 if( stepStatus != fGeomBoundary ) {
714 // NO single scattering in case of skin or << 584 presafety = ComputeSafety(sp->GetPosition(),fTheTrueStepLenght);
715 // model will be single or even no scatter <<
716 // compared to the elastic mean free path. <<
717 // <<
718 // indicate that MSC needs to be done (alw <<
719 fIsMultipleSacettring = true; <<
720 if (stepStatus!=fGeomBoundary) { <<
721 presafety = ComputeSafety(sp->GetPositio <<
722 } <<
723 // Far from boundary-> in optimized mode d <<
724 if ((distance<presafety) && (gIsOptimizati <<
725 fIsNoDisplace = true; <<
726 } else { <<
727 // Urban like <<
728 if (firstStep || (stepStatus==fGeomBound <<
729 rangeinit = currentRange; <<
730 fr = facrange; <<
731 // We don't use this: we won't converge to the <<
732 // decreasing range-factor. <<
733 // rangeinit = std::max(rangeinit <<
734 // if(fLambda1 > lambdalimit) { <<
735 // fr *= (0.75+0.25*fLambda1/la <<
736 // } <<
737 <<
738 } 585 }
739 //step limit << 586
740 tlimit = std::max(fr*rangeinit, facsafet << 587 // Far from boundary-> in optimized mode do not sample dispalcement.
741 // first step randomization << 588 if( (currentRange < presafety) && (fgIsOptimizationOn)) {
742 if (firstStep || stepStatus==fGeomBounda << 589 // fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*mlogAlmostZero);
743 fTheTrueStepLenght = std::min(fTheTrue << 590 // fTheZPathLenght = fLambda1*(1.-G4Exp(-fTheTrueStepLenght/fLambda1));
>> 591 fIsNoDisplace = true;
744 } else { 592 } else {
745 fTheTrueStepLenght = std::min(fTheTrue << 593 // Urban like
>> 594 if(firstStep || (stepStatus == fGeomBoundary)) {
>> 595 rangeinit = currentRange;
>> 596 fr = facrange;
>> 597 rangeinit = std::max(rangeinit, fLambda1);
>> 598 if(fLambda1 > lambdalimit) {
>> 599 fr *= (0.75+0.25*fLambda1/lambdalimit);
>> 600 }
>> 601 }
>> 602
>> 603 //step limit
>> 604 tlimit = std::max(fr*rangeinit, facsafety*presafety);
>> 605
>> 606 // first step randomization
>> 607 if(firstStep || stepStatus == fGeomBoundary) {
>> 608 fTheTrueStepLenght = std::min(fTheTrueStepLenght, Randomizetlimit());
>> 609 } else {
>> 610 fTheTrueStepLenght = std::min(fTheTrueStepLenght, tlimit);
>> 611 }
>> 612
>> 613 // Do not let the true step length to be longer than the one that has a
>> 614 // corresponding geometrical length almost one first transport mean free
>> 615 // path. It is important for the g->t conversion.
>> 616 fTheTrueStepLenght = std::min(fTheTrueStepLenght, fLambda1*mlogAlmostZero);
746 } 617 }
747 } << 618
>> 619 // indicate that MSC needs to be done (always after transportation)
>> 620 fIsMultipleSacettring = true;
748 } 621 }
749 // << 622
750 // unset first-step 623 // unset first-step
751 firstStep =false; 624 firstStep =false;
>> 625
752 // performe single scattering, multiple scat 626 // performe single scattering, multiple scattering if this later can be done safely here
753 if (fIsEverythingWasDone) { << 627 if(fIsEverythingWasDone){
754 if (fIsSingleScattering) { << 628 if(fIsSingleScattering){
755 // sample single scattering 629 // sample single scattering
756 //G4double ekin = 0.5*(currentKinEnerg << 630 G4double cost,sint,cosPhi,sinPhi;
757 G4double lekin = G4Log(currentKinEnergy << 631 SingleScattering(cost, sint);
758 G4double pt2 = currentKinEnergy*(curr << 632 G4double phi = CLHEP::twopi*G4UniformRand();
759 G4double beta2 = pt2/(pt2+CLHEP::electr << 633 sinPhi = std::sin(phi);
760 G4double cost = fGSTable->SingleScatte << 634 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 635 fTheNewDirection.set(sint*cosPhi,sint*sinPhi,cost);
771 } else if (fIsMultipleSacettring) { << 636 } else if(fIsMultipleSacettring){
772 // sample multiple scattering 637 // sample multiple scattering
773 SampleMSC(); // fTheZPathLenght, fTheDis 638 SampleMSC(); // fTheZPathLenght, fTheDisplacementVector and fTheNewDirection will be set
774 } // and if single scattering but it was l 639 } // and if single scattering but it was longer => nothing to do
775 } //else { do nothing here but after transpo 640 } //else { do nothing here but after transportation
776 // << 641
777 return ConvertTrueToGeom(fTheTrueStepLenght, 642 return ConvertTrueToGeom(fTheTrueStepLenght,currentMinimalStep);
778 } 643 }
779 644
780 645
781 G4double G4GoudsmitSaundersonMscModel::Compute << 646 ////////////////////////////////////////////////////////////////////////////////
>> 647 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
>> 648 {
>> 649 //
782 // convert true ->geom 650 // convert true ->geom
783 // It is called from the step limitation Com 651 // It is called from the step limitation ComputeTruePathLengthLimit if
784 // !fIsEverythingWasDone but protect: 652 // !fIsEverythingWasDone but protect:
785 par1 = -1.; << 653 //
786 par2 = par3 = 0.; << 654 par1 = -1. ;
>> 655 par2 = par3 = 0. ;
>> 656
787 // if fIsEverythingWasDone = TRUE => fTheZP 657 // if fIsEverythingWasDone = TRUE => fTheZPathLenght is already set
788 // so return with the already known value 658 // so return with the already known value
789 // Otherwise: 659 // Otherwise:
790 if (!fIsEverythingWasDone) { << 660 if(!fIsEverythingWasDone) {
791 // this correction needed to run MSC with << 661 // this correction needed to run MSC with eIoni and eBrem inactivated
792 // and makes no harm for a normal run << 662 // and makes no harm for a normal run
793 fTheTrueStepLenght = std::min(fTheTrueStep << 663 fTheTrueStepLenght = std::min(fTheTrueStepLenght,currentRange);
794 // do the true -> geom transformation << 664
795 fTheZPathLenght = fTheTrueStepLenght; << 665 // do the true -> geom transformation
796 // z = t for very small true-path-length << 666 fTheZPathLenght = fTheTrueStepLenght;
797 if (fTheTrueStepLenght<tlimitminfix2) { << 667
798 return fTheZPathLenght; << 668 // z = t for very small true-path-length
799 } << 669 if(fTheTrueStepLenght < tlimitminfix2) return fTheZPathLenght;
800 G4double tau = fTheTrueStepLenght/fLambda1 << 670
801 if (tau<=tausmall) { << 671 G4double tau = fTheTrueStepLenght/fLambda1;
802 fTheZPathLenght = std::min(fTheTrueStepL << 672
803 } else if (fTheTrueStepLenght<currentRang << 673 if (tau <= tausmall) {
804 if (tau<taulim) fTheZPathLenght = fTheTr << 674 fTheZPathLenght = std::min(fTheTrueStepLenght, fLambda1);
805 else fTheZPathLenght = fLambd << 675 } else if (fTheTrueStepLenght < currentRange*dtrl) {
806 } else if (currentKinEnergy<mass || fTheTr << 676 if(tau < taulim) fTheZPathLenght = fTheTrueStepLenght*(1.-0.5*tau) ;
807 par1 = 1./currentRange ; // alpha =1 << 677 else fTheZPathLenght = fLambda1*(1.-G4Exp(-tau));
808 par2 = 1./(par1*fLambda1) ; // 1/(alpha << 678 } else if(currentKinEnergy < mass || fTheTrueStepLenght == currentRange) {
809 par3 = 1.+par2 ; // 1+1/ << 679 par1 = 1./currentRange ; // alpha =1/range_init for Ekin<mass
810 if (fTheTrueStepLenght<currentRange) { << 680 par2 = 1./(par1*fLambda1) ; // 1/(alphaxlambda01)
811 fTheZPathLenght = 1./(par1*par3) * (1. << 681 par3 = 1.+par2 ; // 1+1/
812 } else { << 682 if(fTheTrueStepLenght < currentRange) {
813 fTheZPathLenght = 1./(par1*par3); << 683 fTheZPathLenght = 1./(par1*par3) * (1.-std::pow(1.-par1*fTheTrueStepLenght,par3));
814 } << 684 } else {
>> 685 fTheZPathLenght = 1./(par1*par3);
>> 686 }
>> 687
815 } else { 688 } else {
816 G4double rfin = std::max(currentRange << 689 G4double rfin = std::max(currentRange-fTheTrueStepLenght, 0.01*currentRange);
817 G4double T1 = GetEnergy(particle,rf << 690 G4double T1 = GetEnergy(particle,rfin,currentCouple);
818 G4double lambda1 = GetTransportMeanFreeP 691 G4double lambda1 = GetTransportMeanFreePathOnly(particle,T1);
819 // << 692
820 par1 = (fLambda1-lambda1)/(fLambda1*fThe 693 par1 = (fLambda1-lambda1)/(fLambda1*fTheTrueStepLenght); // alpha
821 par2 = 1./(par1*fLambda1); 694 par2 = 1./(par1*fLambda1);
822 par3 = 1.+par2 ; 695 par3 = 1.+par2 ;
823 G4Pow *g4calc = G4Pow::GetInstance(); << 696 G4Pow *g4pow = G4Pow::GetInstance();
824 fTheZPathLenght = 1./(par1*par3) * (1.-g << 697 fTheZPathLenght = 1./(par1*par3) * (1.-g4pow->powA(1.-par1*fTheTrueStepLenght,par3));
825 } 698 }
>> 699
>> 700 fTheZPathLenght = std::min(fTheZPathLenght, fLambda1);
826 } 701 }
827 fTheZPathLenght = std::min(fTheZPathLenght, <<
828 // <<
829 return fTheZPathLenght; <<
830 } <<
831 702
>> 703 return fTheZPathLenght;
>> 704 }
832 705
833 G4double G4GoudsmitSaundersonMscModel::Compute << 706 ////////////////////////////////////////////////////////////////////////////////
>> 707 G4double G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
>> 708 {
834 // init 709 // init
835 fIsEndedUpOnBoundary = false; 710 fIsEndedUpOnBoundary = false;
>> 711
836 // step defined other than transportation 712 // step defined other than transportation
837 if (geomStepLength==fTheZPathLenght) { << 713 if(geomStepLength == fTheZPathLenght) {
838 return fTheTrueStepLenght; 714 return fTheTrueStepLenght;
839 } 715 }
>> 716
840 // else :: 717 // else ::
841 // - set the flag that transportation was th 718 // - set the flag that transportation was the winer so DoNothin in DOIT !!
842 // - convert geom -> true by using the mean 719 // - convert geom -> true by using the mean value
843 fIsEndedUpOnBoundary = true; // OR LAST STEP 720 fIsEndedUpOnBoundary = true; // OR LAST STEP
844 fTheZPathLenght = geomStepLength; << 721
845 // was a short single scattering step << 722 fTheZPathLenght = geomStepLength;
846 if (fIsEverythingWasDone && !fIsMultipleSace << 723
847 fTheTrueStepLenght = geomStepLength; <<
848 return fTheTrueStepLenght; <<
849 } <<
850 // t = z for very small step 724 // t = z for very small step
851 if (geomStepLength<tlimitminfix2) { << 725 if(geomStepLength < tlimitminfix2) {
852 fTheTrueStepLenght = geomStepLength; 726 fTheTrueStepLenght = geomStepLength;
853 // recalculation 727 // recalculation
854 } else { 728 } else {
855 G4double tlength = geomStepLength; 729 G4double tlength = geomStepLength;
856 if (geomStepLength>fLambda1*tausmall) { << 730 if(geomStepLength > fLambda1*tausmall ) {
857 if (par1< 0.) { << 731 if(par1 < 0.) {
858 tlength = -fLambda1*G4Log(1.-geomStepL 732 tlength = -fLambda1*G4Log(1.-geomStepLength/fLambda1) ;
859 } else { 733 } else {
860 if (par1*par3*geomStepLength<1.) { << 734 if(par1*par3*geomStepLength < 1.) {
861 G4Pow *g4calc = G4Pow::GetInstance() << 735 G4Pow *g4pow = G4Pow::GetInstance();
862 tlength = (1.-g4calc->powA( 1.-par1* << 736 tlength = (1.-g4pow->powA( 1.-par1*par3*geomStepLength,1./par3))/par1;
863 } else { 737 } else {
864 tlength = currentRange; 738 tlength = currentRange;
865 } 739 }
866 } 740 }
867 if (tlength<geomStepLength || tlength>fT << 741 if(tlength < geomStepLength) { tlength = geomStepLength; }
868 tlength = geomStepLength; << 742 else if(tlength > fTheTrueStepLenght) { tlength = geomStepLength; }
869 } <<
870 } 743 }
871 fTheTrueStepLenght = tlength; 744 fTheTrueStepLenght = tlength;
872 } 745 }
873 // <<
874 return fTheTrueStepLenght; 746 return fTheTrueStepLenght;
875 } 747 }
876 748
>> 749 ////////////////////////////////////////////////////////////////////////////////
877 G4ThreeVector& 750 G4ThreeVector&
878 G4GoudsmitSaundersonMscModel::SampleScattering << 751 G4GoudsmitSaundersonMscModel::SampleScattering(const G4ThreeVector& oldDirection,
879 if (steppingAlgorithm==fUseDistanceToBoundar << 752 G4double)
880 // single scattering was and scattering ha << 753 {
881 fTheNewDirection.rotateUz(oldDirection); << 754
882 fParticleChange->ProposeMomentumDirection( << 755
883 return fTheDisplacementVector; << 756 if(fIsEndedUpOnBoundary && ( (steppingAlgorithm == fUseSafety) || (steppingAlgorithm == fUseDistanceToBoundary) )){
884 } else if (steppingAlgorithm==fUseSafetyPlus << 757 // transportation was the winer => do nothing
885 if (fIsEndedUpOnBoundary) { // do nothing << 758 fTheDisplacementVector.set(0.,0.,0.);
886 return fTheDisplacementVector; << 759 return fTheDisplacementVector;
887 } else if (fIsEverythingWasDone) { // evry << 760 } else if(fIsEverythingWasDone) {
888 // check single scattering and see if it << 761 // everything was done in the step limit => rotate and return
889 if (fIsSingleScattering) { << 762 // but do it only something happend in the MSC
890 fTheNewDirection.rotateUz(oldDirection << 763 if(!fIsNoScatteringInMSC) {
891 fParticleChange->ProposeMomentumDirect << 764 fTheNewDirection.rotateUz(oldDirection);
892 return fTheDisplacementVector; <<
893 } <<
894 // check if multiple scattering happened <<
895 if (fIsMultipleSacettring && !fIsNoScatt <<
896 fTheNewDirection.rotateUz(oldDirect <<
897 fTheDisplacementVector.rotateUz(old <<
898 fParticleChange->ProposeMomentumDir <<
899 } <<
900 // The only thing that could happen if w <<
901 // is that single scattering was tried <<
902 // So no displacement and no scattering <<
903 return fTheDisplacementVector; <<
904 } <<
905 // <<
906 // The only thing that could still happen <<
907 // optimization branch: so sample MSC angl <<
908 } <<
909 //else MSC needs to be done here <<
910 SampleMSC(); <<
911 if (!fIsNoScatteringInMSC) { <<
912 fTheNewDirection.rotateUz(oldDirection); <<
913 fParticleChange->ProposeMomentumDirection( <<
914 if (!fIsNoDisplace) { <<
915 fTheDisplacementVector.rotateUz(oldDirec 765 fTheDisplacementVector.rotateUz(oldDirection);
916 } << 766 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
917 } << 767 }
918 // << 768 return fTheDisplacementVector;
919 return fTheDisplacementVector; << 769 }
>> 770
>> 771 //else {
>> 772 // MSC needs to be done here: cause we went further than safety but did not hit boundary
>> 773 SampleMSC();
>> 774 if(!fIsNoScatteringInMSC) {
>> 775 fTheNewDirection.rotateUz(oldDirection);
>> 776 fTheDisplacementVector.rotateUz(oldDirection);
>> 777 fParticleChange->ProposeMomentumDirection(fTheNewDirection);
>> 778 }
>> 779 return fTheDisplacementVector;
920 } 780 }
921 781
922 782
923 void G4GoudsmitSaundersonMscModel::SampleMSC() << 783 void G4GoudsmitSaundersonMscModel::SingleScattering(G4double &cost, G4double &sint){
>> 784 G4double rand1 = G4UniformRand();
>> 785 // sampling 1-cos(theta)
>> 786 G4double dum0 = 2.0*fScrA*rand1/(1.0 - rand1 + fScrA);
>> 787 // add protections
>> 788 if(dum0 < 0.0) dum0 = 0.0;
>> 789 else if(dum0 > 2.0 ) dum0 = 2.0;
>> 790 // compute cos(theta) and sin(theta)
>> 791 cost = 1.0 - dum0;
>> 792 sint = std::sqrt(dum0*(2.0 - dum0));
>> 793 }
>> 794
>> 795
>> 796
>> 797 ////////////////////////////////////////////////////////////////////////////////
>> 798 void G4GoudsmitSaundersonMscModel::SampleMSC(){
924 fIsNoScatteringInMSC = false; 799 fIsNoScatteringInMSC = false;
925 // kinetic energy is assumed to be in Geant4 800 // kinetic energy is assumed to be in Geant4 internal energy unit which is MeV
926 G4double kineticEnergy = currentKinEnergy; 801 G4double kineticEnergy = currentKinEnergy;
927 // << 802
>> 803 ///////////////////////////////////////////
928 // Energy loss correction: 2 version 804 // Energy loss correction: 2 version
929 G4double eloss = 0.0; 805 G4double eloss = 0.0;
930 // if (fTheTrueStepLenght > currentRange*dtrl << 806 if (fTheTrueStepLenght > currentRange*dtrl) {
931 eloss = kineticEnergy - GetEnergy(particle,c << 807 eloss = kineticEnergy -
932 // } else { << 808 GetEnergy(particle,currentRange-fTheTrueStepLenght,currentCouple);
933 // eloss = fTheTrueStepLenght*GetDEDX(parti << 809 } else {
934 // } << 810 eloss = fTheTrueStepLenght*GetDEDX(particle,kineticEnergy,currentCouple);
>> 811 }
>> 812
935 813
936 G4double tau = 0.;// = kineticEnergy/ele 814 G4double tau = 0.;// = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
937 G4double tau2 = 0.;// = tau*tau; 815 G4double tau2 = 0.;// = tau*tau;
938 G4double eps0 = 0.;// = eloss/kineticEnerg 816 G4double eps0 = 0.;// = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
939 G4double epsm = 0.;// = eloss/kineticEnerg 817 G4double epsm = 0.;// = eloss/kineticEnergy; // energy loss fraction to the mean step energy
940 818
>> 819
941 // - init. 820 // - init.
942 G4double efEnergy = kineticEnergy; 821 G4double efEnergy = kineticEnergy;
943 G4double efStep = fTheTrueStepLenght; 822 G4double efStep = fTheTrueStepLenght;
944 823
945 G4double kineticEnergy0 = kineticEnergy; << 824 if(fgIsUseAccurate) { // - use accurate energy loss correction
946 if (gIsUseAccurate) { // - use accurate e << 825 G4double kineticEnergy0 = kineticEnergy;
947 kineticEnergy -= 0.5*eloss; // mean ener << 826 kineticEnergy -= 0.5*eloss; // mean energy along the full step
948 // other parameters for energy loss correc << 827
949 tau = kineticEnergy/electron_m << 828 // other parameters for energy loss corrections
950 tau2 = tau*tau; << 829 tau = kineticEnergy/electron_mass_c2; // where kinEnergy is the mean kinetic energy
951 eps0 = eloss/kineticEnergy0; // << 830 tau2 = tau*tau;
952 epsm = eloss/kineticEnergy; // << 831 eps0 = eloss/kineticEnergy0; // energy loss fraction to the begin step energy
953 << 832 epsm = eloss/kineticEnergy; // energy loss fraction to the mean step energy
954 efEnergy = kineticEnergy * (1.-epsm << 833 // efEnergy = kineticEnergy0 * (1.0 - epsm*(6.0+10.0*tau+5.0*tau2)/(24.0*tau2+72.0*tau+48.0));
955 G4double dum = 0.166666*(4.+tau*(6.+tau << 834 efEnergy = kineticEnergy0 * ( 1.0 - 0.5*eps0 - eps0*eps0/(12.0*(2.0-eps0))* ( (5.*tau2 +10.*tau +6.)/( (tau+1.)*(tau+2.) ) ) );
956 efStep = fTheTrueStepLenght*(1.-d << 835
>> 836 // the real efStep will be s*(1-efStep)
>> 837 // 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);
>> 838
>> 839 // efStep = 0.166666*eps0*epsm*(4.0+tau*(6.0+tau*(7.0+tau*(4.0+tau))));
>> 840 // efStep /= ((tau2+3.0*tau+2)*(tau2+3.0*tau+2));
>> 841 // efStep = fTheTrueStepLenght*(1.0-efStep);
>> 842 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.)) ) );
>> 843 efStep *=fTheTrueStepLenght;
957 } else { // - t 844 } else { // - take only mean energy
958 kineticEnergy -= 0.5*eloss; // mean ener << 845 kineticEnergy -= 0.5*eloss; // mean energy along the full step
959 efEnergy = kineticEnergy; << 846 efEnergy = kineticEnergy;
960 G4double factor = 1./(1.+0.9784671*kinetic << 847 efStep = fTheTrueStepLenght;
961 eps0 = eloss/kineticEnergy0; <<
962 epsm = eps0/(1.-0.5*eps0); <<
963 G4double temp = 0.3*(1 -factor*(1.-0.333 <<
964 efStep = fTheTrueStepLenght*(1.+t <<
965 } 848 }
966 // << 849 /////////////////////////////////////////////////////////////////////////
967 // compute elastic mfp, first transport mfp, << 850
968 // if it was requested by the user) << 851
969 fLambda1 = GetTransportMeanFreePath(particle << 852 ///////////////////////////////////////////
970 // s/lambda_el << 853 // Compute elastic mfp, first transport mfp, screening parameter, and G1
>> 854 fLambda1 = GetTransportMeanFreePath(particle,efEnergy);
>> 855
971 G4double lambdan=0.; 856 G4double lambdan=0.;
972 if (fLambda0>0.0) { << 857
973 lambdan=efStep/fLambda0; << 858 if(fLambda0>0.0) { lambdan=efStep/fLambda0; }
974 } << 859 if(lambdan<=1.0e-12) {
975 if (lambdan<=1.0e-12) { << 860 if(fIsEverythingWasDone)
976 if (fIsEverythingWasDone) { <<
977 fTheZPathLenght = fTheTrueStepLenght; 861 fTheZPathLenght = fTheTrueStepLenght;
978 } <<
979 fIsNoScatteringInMSC = true; 862 fIsNoScatteringInMSC = true;
980 return; 863 return;
981 } 864 }
982 // first moment: 2.* lambdan *scrA*((1.+scrA << 865
983 G4double Qn1 = lambdan *fG1; << 866 G4double Qn1 = lambdan *fG1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
984 // sample scattering angles << 867
985 // new direction, relative to the orriginal << 868 /////////////////////////////////////////////////////////////////////////////
986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, c << 869 // Sample scattering angles
987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, c << 870 // new direction, relative to the orriginal one is in {us,vs,ws}
988 G4double uss = 0.0, vss = 0.0, w << 871 G4double cosTheta1 =1.0,sinTheta1 =0.0, cosTheta2 =1.0, sinTheta2 =0.0;
989 G4double x_coord = 0.0, y_coord = 0.0, z << 872 G4double cosPhi1=1.0,sinPhi1=0.0, cosPhi2 =1.0, sinPhi2 =0.0;
990 G4double u2 = 0.0, v2 = 0.0; << 873 G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
>> 874 G4double u2=0.0,v2=0.0;
>> 875
991 // if we are above the upper grid limit with 876 // if we are above the upper grid limit with lambdaxG1=true-length/first-trans-mfp
992 // => izotropic distribution: lambG1_max =7. << 877 // => izotropic distribution: lambG1_max =7.992
993 if (0.5*Qn1 > 7.0){ << 878 if(0.5*Qn1 > 7.99){
994 cosTheta1 = 1.-2.*G4UniformRand(); << 879 cosTheta1 = 1.-2.*G4UniformRand();
995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+c << 880 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+cosTheta1));
996 cosTheta2 = 1.-2.*G4UniformRand(); << 881 G4double phi1 = CLHEP::twopi*G4UniformRand();
997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+c << 882 sinPhi1 = std::sin(phi1);
>> 883 cosPhi1 = std::cos(phi1);
>> 884 us = sinTheta1*cosPhi1;
>> 885 vs = sinTheta1*sinPhi1;
>> 886 ws = cosTheta1;
998 } else { 887 } else {
999 // sample 2 scattering cost1, sint1, cost 888 // sample 2 scattering cost1, sint1, cost2 and sint2 for half path
1000 G4double lekin = G4Log(efEnergy); << 889 fgGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta1, sinTheta1);
1001 G4double pt2 = efEnergy*(efEnergy+2.0 << 890 fgGSTable->Sampling(0.5*lambdan, 0.5*Qn1, fScrA, cosTheta2, sinTheta2);
1002 G4double beta2 = pt2/(pt2+CLHEP::electr << 891 if(cosTheta1 + cosTheta2 == 2.) { // no scattering happened
1003 // backup GS angular dtr pointer (kineti << 892 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 893 fTheZPathLenght = fTheTrueStepLenght;
1017 fIsNoScatteringInMSC = true; 894 fIsNoScatteringInMSC = true;
1018 return; 895 return;
1019 } 896 }
>> 897
>> 898 // sample 2 azimuthal angles
>> 899 G4double phi1 = CLHEP::twopi*G4UniformRand();
>> 900 sinPhi1 = std::sin(phi1);
>> 901 cosPhi1 = std::cos(phi1);
>> 902 G4double phi2 = CLHEP::twopi*G4UniformRand();
>> 903 sinPhi2 = std::sin(phi2);
>> 904 cosPhi2 = std::cos(phi2);
>> 905
>> 906 // compute final direction realtive to z-dir
>> 907 u2 = sinTheta2*cosPhi2;
>> 908 v2 = sinTheta2*sinPhi2;
>> 909 G4double u2p = cosTheta1*u2 + sinTheta1*cosTheta2;
>> 910 us = u2p*cosPhi1 - v2*sinPhi1;
>> 911 vs = u2p*sinPhi1 + v2*cosPhi1;
>> 912 ws = cosTheta1*cosTheta2 - sinTheta1*u2;
1020 } 913 }
1021 // sample 2 azimuthal angles <<
1022 G4double phi1 = CLHEP::twopi*G4UniformRand( <<
1023 sinPhi1 = std::sin(phi1); <<
1024 cosPhi1 = std::cos(phi1); <<
1025 G4double phi2 = CLHEP::twopi*G4UniformRand( <<
1026 sinPhi2 = std::sin(phi2); <<
1027 cosPhi2 = std::cos(phi2); <<
1028 <<
1029 // compute final direction realtive to z-di <<
1030 u2 = sinTheta2*cosPhi2; <<
1031 v2 = sinTheta2*sinPhi2; <<
1032 G4double u2p = cosTheta1*u2 + sinTheta1*cos <<
1033 uss = u2p*cosPhi1 - v2*sinPhi1; <<
1034 vss = u2p*sinPhi1 + v2*cosPhi1; <<
1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2; <<
1036 914
1037 // set new direction (is scattering frame) << 915 ////////////////////////////////////////////////////////////////////
1038 fTheNewDirection.set(uss,vss,wss); << 916 fTheNewDirection.set(us,vs,ws);
1039 917
1040 // set the fTheZPathLenght if we don't samp 918 // set the fTheZPathLenght if we don't sample displacement and
1041 // we should do everything at the step-limi 919 // we should do everything at the step-limit-phase before we return
1042 if(fIsNoDisplace && fIsEverythingWasDone) 920 if(fIsNoDisplace && fIsEverythingWasDone)
1043 fTheZPathLenght = fTheTrueStepLenght; 921 fTheZPathLenght = fTheTrueStepLenght;
1044 922
1045 // in optimized-mode if the current-safety 923 // in optimized-mode if the current-safety > current-range we do not use dispalcement
1046 if(fIsNoDisplace) 924 if(fIsNoDisplace)
1047 return; 925 return;
1048 926
1049 /////////////////////////////////////////// 927 //////////////////////////////////////////////////////////////////////
1050 // Compute final position 928 // Compute final position
1051 Qn1 *= fMCtoQ1; << 929 if(fgIsUseAccurate) {
1052 if (gIsUseAccurate) { << 930 // correction parameters
1053 // correction parameter << 931
1054 G4double par =1.; << 932 G4double par =1.;
1055 if(Qn1<0.7) par = 1.; << 933 if(Qn1<1.0) par = 1.;
1056 else if (Qn1<7.0) par = -0.031376*Qn1+1. << 934 else if (Qn1<7.8) par = -0.031376*Qn1+1.01356;
1057 else par = 0.79; << 935 else par = 0.76;
>> 936
>> 937 //
>> 938 // Compute elastic mfp, first transport mfp, screening parameter, and G1
>> 939 // for the initial energy
>> 940 fLambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
1058 941
1059 // Moments with energy loss correction 942 // Moments with energy loss correction
1060 // --first the uncorrected (for energy l 943 // --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 944 // 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); 945 G4double loga = G4Log(1.0+1.0/fScrA);
1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA 946 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 947 // 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 << 948 G4double rand1 = G4UniformRand();
>> 949 G4double eta = std::sqrt(rand1);
1067 G4double eta1 = 0.5*(1 - eta); // use 950 G4double eta1 = 0.5*(1 - eta); // used more than once
1068 // 0.5 +sqrt(6)/6 = 0.9082483; 951 // 0.5 +sqrt(6)/6 = 0.9082483;
1069 // 1/(4*sqrt(6)) = 0.1020621; 952 // 1/(4*sqrt(6)) = 0.1020621;
1070 // (4-sqrt(6)/(24*sqrt(6))) = 0.02637471 953 // (4-sqrt(6)/(24*sqrt(6))) = 0.026374715
1071 // delta = 0.9082483-(0.1020621-0.026374 954 // delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1 without energy loss cor.
>> 955 Qn1=fTheTrueStepLenght/fLambda0*fG1; // without energy loss
1072 G4double delta = 0.9082483-(0.1020621-0 956 G4double delta = 0.9082483-(0.1020621-0.0263747*gamma)*Qn1;
1073 957
1074 // compute alpha1 and alpha2 for energy 958 // compute alpha1 and alpha2 for energy loss correction
1075 G4double temp1 = 2.0 + tau; 959 G4double temp1 = 2.0 + tau;
1076 G4double temp = (2.0+tau*temp1)/((tau+1 960 G4double temp = (2.0+tau*temp1)/((tau+1.0)*temp1);
1077 //Take logarithmic dependence 961 //Take logarithmic dependence
1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga 962 temp = temp - (tau+1.0)/((tau+2.0)*(loga*(1.0+fScrA)-1.0));
1079 temp = temp * epsm; 963 temp = temp * epsm;
1080 temp1 = 1.0 - temp; 964 temp1 = 1.0 - temp;
1081 delta = delta + 0.40824829*(eps0*(tau+1. 965 delta = delta + 0.40824829*(eps0*(tau+1.0)/((tau+2.0)*
1082 (loga*(1.0+fScrA)-1.0)*(loga*(1. 966 (loga*(1.0+fScrA)-1.0)*(loga*(1.0+2.0*fScrA)-2.0)) - 0.25*temp*temp);
1083 G4double b = eta*delta; 967 G4double b = eta*delta;
1084 G4double c = eta*(1.0-delta); 968 G4double c = eta*(1.0-delta);
1085 969
1086 //Calculate transport direction cosines: 970 //Calculate transport direction cosines:
1087 // ut,vt,wt is the final position divide 971 // ut,vt,wt is the final position divided by the true step length
1088 G4double w1v2 = cosTheta1*v2; 972 G4double w1v2 = cosTheta1*v2;
1089 G4double ut = b*sinTheta1*cosPhi1 + c* << 973 G4double ut = b*sinTheta1*cosPhi1 + c*(cosPhi1*u2 - sinPhi1*w1v2) + eta1*us*temp1;
1090 G4double vt = b*sinTheta1*sinPhi1 + c* << 974 G4double vt = b*sinTheta1*sinPhi1 + c*(sinPhi1*u2 + cosPhi1*w1v2) + eta1*vs*temp1;
1091 G4double wt = eta1*(1+temp) + b* << 975 G4double wt = eta1*(1+temp) + b*cosTheta1 + c*cosTheta2 + eta1*ws*temp1;
1092 976
1093 // long step correction 977 // long step correction
1094 ut *=par; 978 ut *=par;
1095 vt *=par; 979 vt *=par;
1096 wt *=par; 980 wt *=par;
1097 981
1098 // final position relative to the pre-st 982 // final position relative to the pre-step point in the scattering frame
1099 // ut = x_f/s so needs to multiply by s 983 // ut = x_f/s so needs to multiply by s
1100 x_coord = ut*fTheTrueStepLenght; 984 x_coord = ut*fTheTrueStepLenght;
1101 y_coord = vt*fTheTrueStepLenght; 985 y_coord = vt*fTheTrueStepLenght;
1102 z_coord = wt*fTheTrueStepLenght; 986 z_coord = wt*fTheTrueStepLenght;
1103 987
1104 if(fIsEverythingWasDone){ 988 if(fIsEverythingWasDone){
1105 // We sample in the step limit so set 989 // We sample in the step limit so set fTheZPathLenght = transportDistance
1106 // and lateral displacement (x_coord,y 990 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1107 //Calculate transport distance 991 //Calculate transport distance
1108 G4double transportDistance = std::sqr 992 G4double transportDistance = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
1109 // protection 993 // protection
1110 if(transportDistance>fTheTrueStepLengh 994 if(transportDistance>fTheTrueStepLenght)
1111 transportDistance = fTheTrueStepLen 995 transportDistance = fTheTrueStepLenght;
1112 fTheZPathLenght = transportDistance; 996 fTheZPathLenght = transportDistance;
1113 } 997 }
1114 // else:: we sample in the DoIt so 998 // else:: we sample in the DoIt so
1115 // the fTheZPathLenght was already 999 // the fTheZPathLenght was already set and was taken as transport along zet
1116 fTheDisplacementVector.set(x_coord,y_coo 1000 fTheDisplacementVector.set(x_coord,y_coord,z_coord-fTheZPathLenght);
1117 } else { 1001 } else {
1118 // compute zz = <z>/tPathLength 1002 // compute zz = <z>/tPathLength
1119 // s -> true-path-length 1003 // s -> true-path-length
1120 // z -> geom-path-length:: when PRESTA i 1004 // z -> geom-path-length:: when PRESTA is used z =(def.) <z>
1121 // r -> lateral displacement = s/2 sin(t 1005 // r -> lateral displacement = s/2 sin(theta) => x_f = r cos(phi); y_f = r sin(phi)
1122 G4double zz = 0.0; 1006 G4double zz = 0.0;
1123 if(fIsEverythingWasDone){ 1007 if(fIsEverythingWasDone){
1124 // We sample in the step limit so set 1008 // We sample in the step limit so set fTheZPathLenght = transportDistance
1125 // and lateral displacement (x_coord, 1009 // and lateral displacement (x_coord,y_coord,z_coord-transportDistance)
1126 if(Qn1<0.1) { // use 3-order Taylor a 1010 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 1011 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666667 - 0.041666667*Qn1)); // 1/6 =0.166..7 ; 1/24=0.041..
1128 } else { 1012 } else {
1129 zz = (1.-G4Exp(-Qn1))/Qn1; 1013 zz = (1.-G4Exp(-Qn1))/Qn1;
1130 } 1014 }
1131 } else { 1015 } else {
1132 // we sample in the DoIt so 1016 // we sample in the DoIt so
1133 // the fTheZPathLenght was already se 1017 // the fTheZPathLenght was already set and was taken as transport along zet
1134 zz = fTheZPathLenght/fTheTrueStepLeng 1018 zz = fTheZPathLenght/fTheTrueStepLenght;
1135 } 1019 }
1136 1020
1137 G4double rr = (1.-zz*zz)/(1.-wss*wss); / << 1021 G4double rr = (1.-zz*zz)/(1.-ws*ws); // s^2 >= <z>^2+r^2 :: where r^2 = s^2/4 sin^2(theta)
1138 if(rr >= 0.25) rr = 0.25; // 1022 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: 1023 G4double rperp = fTheTrueStepLenght*std::sqrt(rr); // this is r/sint
1140 x_coord = rperp*uss; << 1024 x_coord = rperp*us;
1141 y_coord = rperp*vss; << 1025 y_coord = rperp*vs;
1142 z_coord = zz*fTheTrueStepLenght; 1026 z_coord = zz*fTheTrueStepLenght;
1143 1027
1144 if(fIsEverythingWasDone){ 1028 if(fIsEverythingWasDone){
1145 G4double transportDistance = std::sqrt 1029 G4double transportDistance = std::sqrt(x_coord*x_coord + y_coord*y_coord + z_coord*z_coord);
1146 fTheZPathLenght = transportDistance; 1030 fTheZPathLenght = transportDistance;
1147 } 1031 }
1148 1032
1149 fTheDisplacementVector.set(x_coord,y_coo 1033 fTheDisplacementVector.set(x_coord,y_coord,z_coord- fTheZPathLenght);
1150 } 1034 }
1151 } 1035 }
1152 1036