Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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