Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.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 ]

Diff markup

Differences between /processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc (Version 11.3.0) and /processes/electromagnetic/standard/src/G4GoudsmitSaundersonTable.cc (Version 11.1.2)


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