| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
>> 26 // $Id: G4GoudsmitSaundersonTable.cc 94933 2015-12-18 09:22:52Z gcosmo $
26 // 27 //
27 // ------------------------------------------- 28 // -----------------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class implementation file 30 // GEANT4 Class implementation file
30 // 31 //
31 // File name: G4GoudsmitSaundersonTable 32 // File name: G4GoudsmitSaundersonTable
32 // 33 //
33 // Author: Mihaly Novak / (Omrane Kadri 34 // Author: Mihaly Novak / (Omrane Kadri)
34 // 35 //
35 // Creation date: 20.02.2009 36 // Creation date: 20.02.2009
36 // 37 //
37 // Class description: 38 // Class description:
38 // Class to handle multiple scattering angul 39 // Class to handle multiple scattering angular distributions precomputed by
39 // using Kawrakow-Bielajew Goudsmit-Saunders 40 // using Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened
40 // Rutherford DCS for elastic scattering of 41 // Rutherford DCS for elastic scattering of electrons/positrons [1,2]. This
41 // class is used by G4GoudsmitSaundersonMscM 42 // class is used by G4GoudsmitSaundersonMscModel to sample the angular
42 // deflection of electrons/positrons after t 43 // deflection of electrons/positrons after travelling a given path.
43 // 44 //
44 // Modifications: 45 // Modifications:
45 // 04.03.2009 V.Ivanchenko cleanup and format 46 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
46 // 26.08.2009 O.Kadri: avoiding unuseful calcu 47 // 26.08.2009 O.Kadri: avoiding unuseful calculations and optimizing the root
47 // finding parameter error 48 // finding parameter error's within SampleTheta method
48 // 08.02.2010 O.Kadri: reduce delared variable 49 // 08.02.2010 O.Kadri: reduce delared variables; reduce error of finding root
49 // in secant method 50 // in secant method
50 // 26.03.2010 O.Kadri: minimum of used arrays 51 // 26.03.2010 O.Kadri: minimum of used arrays in computation within the dichotomie
51 // finding method the erro 52 // finding method the error was the lowest value of uvalues
52 // 12.05.2010 O.Kadri: changing of sqrt((b-a)* 53 // 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 54 // 18.05.2015 M. Novak This class has been completely replaced (only the original
54 // class name was kept; class descr 55 // class name was kept; class description was also inserted):
55 // A new version of Kawrakow-Bielaj 56 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model
56 // based on the screened Rutherford 57 // based on the screened Rutherford DCS for elastic scattering of
57 // electrons/positrons has been int 58 // electrons/positrons has been introduced[1,2]. The corresponding MSC
58 // angular distributions over a 2D 59 // angular distributions over a 2D parameter grid have been recomputed
59 // and the CDFs are now stored in a 60 // and the CDFs are now stored in a variable transformed (smooth) form
60 // together with the corresponding 61 // together with the corresponding rational interpolation parameters.
61 // The new version is several times 62 // The new version is several times faster, more robust and accurate
62 // compared to the earlier version 63 // compared to the earlier version (G4GoudsmitSaundersonMscModel class
63 // that use these data has been als 64 // that use these data has been also completely replaced)
64 // 28.04.2017 M. Novak: New representation of <<
65 // significantly reduced data size. <<
66 // 23.08.2017 M. Novak: Added funtionality to <<
67 // base GS angular distributions an <<
68 // parameter, first and second mome <<
69 // activated in the GS-MSC model. <<
70 // 65 //
71 // References: 66 // References:
72 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-20 67 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-208
73 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19 68 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336
74 // 69 //
75 // ------------------------------------------- 70 // -----------------------------------------------------------------------------
76 71
77 #include "G4GoudsmitSaundersonTable.hh" 72 #include "G4GoudsmitSaundersonTable.hh"
78 73
79 <<
80 #include "G4PhysicalConstants.hh" <<
81 #include "Randomize.hh" 74 #include "Randomize.hh"
82 #include "G4Log.hh" << 75 #include "G4PhysicalConstants.hh"
83 #include "G4Exp.hh" <<
84 <<
85 #include "G4GSMottCorrection.hh" <<
86 #include "G4MaterialTable.hh" 76 #include "G4MaterialTable.hh"
87 #include "G4Material.hh" 77 #include "G4Material.hh"
88 #include "G4MaterialCutsCouple.hh" << 78 #include "G4Log.hh"
89 #include "G4ProductionCutsTable.hh" << 79 #include "G4Exp.hh"
90 #include "G4EmParameters.hh" <<
91 <<
92 #include "G4String.hh" <<
93 80
94 #include <fstream> 81 #include <fstream>
95 #include <cstdlib> 82 #include <cstdlib>
>> 83 #include <cstdio>
96 #include <cmath> 84 #include <cmath>
97 85
98 #include <iostream> <<
99 #include <iomanip> <<
100 86
101 // perecomputed GS angular distributions, base << 87 const G4double G4GoudsmitSaundersonTable::fgLambdaValues[] ={
102 // are the same for e- and e+ so make sure we << 88 1.000000000000e+00, 1.165914401180e+00, 1.359356390879e+00, 1.584893192461e+00,
103 G4bool G4GoudsmitSaundersonTable::gIsInitialis << 89 1.847849797422e+00, 2.154434690032e+00, 2.511886431510e+00, 2.928644564625e+00,
104 // << 90 3.414548873834e+00, 3.981071705535e+00, 4.641588833613e+00, 5.411695265465e+00,
105 std::vector<G4GoudsmitSaundersonTable::GSMSCAn << 91 6.309573444802e+00, 7.356422544596e+00, 8.576958985909e+00, 1.000000000000e+01,
106 std::vector<G4GoudsmitSaundersonTable::GSMSCAn << 92 1.165914401180e+01, 1.359356390879e+01, 1.584893192461e+01, 1.847849797422e+01,
107 // << 93 2.154434690032e+01, 2.511886431510e+01, 2.928644564625e+01, 3.414548873834e+01,
108 std::vector<double> G4GoudsmitSaundersonTable: << 94 3.981071705535e+01, 4.641588833613e+01, 5.411695265465e+01, 6.309573444802e+01,
109 std::vector<double> G4GoudsmitSaundersonTable: << 95 7.356422544596e+01, 8.576958985909e+01, 1.000000000000e+02, 1.165914401180e+02,
>> 96 1.359356390879e+02, 1.584893192461e+02, 1.847849797422e+02, 2.154434690032e+02,
>> 97 2.511886431510e+02, 2.928644564625e+02, 3.414548873834e+02, 3.981071705535e+02,
>> 98 4.641588833613e+02, 5.411695265465e+02, 6.309573444802e+02, 7.356422544596e+02,
>> 99 8.576958985909e+02, 1.000000000000e+03, 1.165914401180e+03, 1.359356390879e+03,
>> 100 1.584893192461e+03, 1.847849797422e+03, 2.154434690032e+03, 2.511886431510e+03,
>> 101 2.928644564625e+03, 3.414548873834e+03, 3.981071705535e+03, 4.641588833613e+03,
>> 102 5.411695265465e+03, 6.309573444802e+03, 7.356422544596e+03, 8.576958985909e+03,
>> 103 1.000000000000e+04, 1.165914401180e+04, 1.359356390879e+04, 1.584893192461e+04,
>> 104 1.847849797422e+04, 2.154434690032e+04, 2.511886431510e+04, 2.928644564625e+04,
>> 105 3.414548873834e+04, 3.981071705535e+04, 4.641588833613e+04, 5.411695265465e+04,
>> 106 6.309573444802e+04, 7.356422544596e+04, 8.576958985909e+04, 1.000000000000e+05
>> 107 };
>> 108
>> 109 const G4double G4GoudsmitSaundersonTable::fgLamG1Values[]={
>> 110 0.0010, 0.0509, 0.1008, 0.1507, 0.2006, 0.2505, 0.3004, 0.3503, 0.4002,
>> 111 0.4501, 0.5000, 0.5499, 0.5998, 0.6497, 0.6996, 0.7495, 0.7994, 0.8493,
>> 112 0.8992, 0.9491, 0.9990
>> 113 };
>> 114 const G4double G4GoudsmitSaundersonTable::fgLamG1ValuesII[]={
>> 115 0.999, 1.332, 1.665, 1.998, 2.331, 2.664, 2.997, 3.33, 3.663,
>> 116 3.996, 4.329, 4.662, 4.995, 5.328, 5.661, 5.994, 6.327,
>> 117 6.66, 6.993, 7.326, 7.659, 7.992
>> 118 };
110 119
>> 120 const G4double G4GoudsmitSaundersonTable::fgUValues[]={
>> 121 0.00, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
>> 122 0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19,
>> 123 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29,
>> 124 0.30, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39,
>> 125 0.40, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49,
>> 126 0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59,
>> 127 0.60, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69,
>> 128 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79,
>> 129 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89,
>> 130 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99,
>> 131 1.00
>> 132 };
>> 133
>> 134 const G4double G4GoudsmitSaundersonTable::fgG1Values[]={
>> 135 1.00000000000000e-10, 1.15478198468946e-10, 1.33352143216332e-10, 1.53992652605949e-10, 1.77827941003892e-10,
>> 136 2.05352502645715e-10, 2.37137370566166e-10, 2.73841963426436e-10, 3.16227766016838e-10, 3.65174127254838e-10,
>> 137 4.21696503428582e-10, 4.86967525165863e-10, 5.62341325190349e-10, 6.49381631576211e-10, 7.49894209332456e-10,
>> 138 8.65964323360065e-10, 1.00000000000000e-09, 1.15478198468946e-09, 1.33352143216332e-09, 1.53992652605949e-09,
>> 139 1.77827941003892e-09, 2.05352502645715e-09, 2.37137370566166e-09, 2.73841963426436e-09, 3.16227766016838e-09,
>> 140 3.65174127254838e-09, 4.21696503428582e-09, 4.86967525165863e-09, 5.62341325190349e-09, 6.49381631576211e-09,
>> 141 7.49894209332456e-09, 8.65964323360065e-09, 1.00000000000000e-08, 1.15478198468946e-08, 1.33352143216332e-08,
>> 142 1.53992652605949e-08, 1.77827941003892e-08, 2.05352502645715e-08, 2.37137370566166e-08, 2.73841963426436e-08,
>> 143 3.16227766016838e-08, 3.65174127254838e-08, 4.21696503428582e-08, 4.86967525165863e-08, 5.62341325190349e-08,
>> 144 6.49381631576211e-08, 7.49894209332456e-08, 8.65964323360065e-08, 1.00000000000000e-07, 1.15478198468946e-07,
>> 145 1.33352143216332e-07, 1.53992652605949e-07, 1.77827941003892e-07, 2.05352502645715e-07, 2.37137370566166e-07,
>> 146 2.73841963426436e-07, 3.16227766016838e-07, 3.65174127254838e-07, 4.21696503428582e-07, 4.86967525165863e-07,
>> 147 5.62341325190349e-07, 6.49381631576211e-07, 7.49894209332456e-07, 8.65964323360065e-07, 1.00000000000000e-06,
>> 148 1.15478198468946e-06, 1.33352143216332e-06, 1.53992652605949e-06, 1.77827941003892e-06, 2.05352502645715e-06,
>> 149 2.37137370566166e-06, 2.73841963426436e-06, 3.16227766016838e-06, 3.65174127254838e-06, 4.21696503428582e-06,
>> 150 4.86967525165863e-06, 5.62341325190349e-06, 6.49381631576211e-06, 7.49894209332456e-06, 8.65964323360065e-06,
>> 151 1.00000000000000e-05, 1.15478198468946e-05, 1.33352143216332e-05, 1.53992652605949e-05, 1.77827941003892e-05,
>> 152 2.05352502645715e-05, 2.37137370566166e-05, 2.73841963426436e-05, 3.16227766016838e-05, 3.65174127254838e-05,
>> 153 4.21696503428582e-05, 4.86967525165863e-05, 5.62341325190349e-05, 6.49381631576211e-05, 7.49894209332456e-05,
>> 154 8.65964323360065e-05, 1.00000000000000e-04, 1.15478198468946e-04, 1.33352143216332e-04, 1.53992652605949e-04,
>> 155 1.77827941003892e-04, 2.05352502645715e-04, 2.37137370566166e-04, 2.73841963426436e-04, 3.16227766016838e-04,
>> 156 3.65174127254838e-04, 4.21696503428582e-04, 4.86967525165863e-04, 5.62341325190349e-04, 6.49381631576211e-04,
>> 157 7.49894209332456e-04, 8.65964323360065e-04, 1.00000000000000e-03, 1.15478198468946e-03, 1.33352143216332e-03,
>> 158 1.53992652605949e-03, 1.77827941003892e-03, 2.05352502645715e-03, 2.37137370566166e-03, 2.73841963426436e-03,
>> 159 3.16227766016838e-03, 3.65174127254838e-03, 4.21696503428582e-03, 4.86967525165863e-03, 5.62341325190349e-03,
>> 160 6.49381631576211e-03, 7.49894209332456e-03, 8.65964323360065e-03, 1.00000000000000e-02, 1.15478198468946e-02,
>> 161 1.33352143216332e-02, 1.53992652605949e-02, 1.77827941003892e-02, 2.05352502645715e-02, 2.37137370566166e-02,
>> 162 2.73841963426436e-02, 3.16227766016838e-02, 3.65174127254838e-02, 4.21696503428582e-02, 4.86967525165863e-02,
>> 163 5.62341325190349e-02, 6.49381631576211e-02, 7.49894209332456e-02, 8.65964323360065e-02, 1.00000000000000e-01,
>> 164 1.15478198468946e-01, 1.33352143216332e-01, 1.53992652605949e-01, 1.77827941003892e-01, 2.05352502645715e-01,
>> 165 2.37137370566166e-01, 2.73841963426436e-01, 3.16227766016838e-01, 3.65174127254838e-01, 4.21696503428582e-01,
>> 166 4.86967525165863e-01, 5.62341325190349e-01, 6.49381631576211e-01, 7.49894209332456e-01, 8.65964323360065e-01
>> 167 };
>> 168
>> 169 const G4double G4GoudsmitSaundersonTable::fgScreeningParam[]={
>> 170 1.92484052135703e-12, 2.23565438533804e-12, 2.59674772669370e-12, 3.01627006395395e-12, 3.50369464193945e-12,
>> 171 4.07003394459337e-12, 4.72809038421179e-12, 5.49274792465712e-12, 6.38131134142232e-12, 7.41390092242269e-12,
>> 172 8.61391169587491e-12, 1.00085477656079e-11, 1.16294440746525e-11, 1.35133899458129e-11, 1.57031711107699e-11,
>> 173 1.82485496926619e-11, 2.12074048158664e-11, 2.46470602564846e-11, 2.86458299060786e-11, 3.32948169025090e-11,
>> 174 3.87000082054842e-11, 4.49847133009623e-11, 5.22924037716594e-11, 6.07900198618435e-11, 7.06718211166503e-11,
>> 175 8.21638709501378e-11, 9.55292598967889e-11, 1.11074189683990e-10, 1.29155060543825e-10, 1.50186727847036e-10,
>> 176 1.74652121757794e-10, 2.03113455838333e-10, 2.36225288153009e-10, 2.74749742338471e-10, 3.19574247380381e-10,
>> 177 3.71732214707185e-10, 4.32427141127727e-10, 5.03060707798248e-10, 5.85265540790081e-10, 6.80943410264647e-10,
>> 178 7.92309775465429e-10, 9.21945734889531e-10, 1.07285861883013e-09, 1.24855266934879e-09, 1.45311149575389e-09,
>> 179 1.69129427781576e-09, 1.96864802125653e-09, 2.29163855873534e-09, 2.66780344424806e-09, 3.10593042087591e-09,
>> 180 3.61626576440636e-09, 4.21075753406168e-09, 4.90333961465076e-09, 5.71026343332307e-09, 6.65048540388811e-09,
>> 181 7.74611952188878e-09, 9.02296613895651e-09, 1.05111298261663e-08, 1.22457414410270e-08, 1.42678020976703e-08,
>> 182 1.66251697709241e-08, 1.93737128201322e-08, 2.25786588894194e-08, 2.63161725354421e-08, 3.06752006784906e-08,
>> 183 3.57596317176864e-08, 4.16908220722131e-08, 4.86105532157964e-08, 5.66844932060409e-08, 6.61062495628185e-08,
>> 184 7.71021154618766e-08, 8.99366289841022e-08, 1.04919086073377e-07, 1.22411172469263e-07, 1.42835908860059e-07,
>> 185 1.66688137634119e-07, 1.94546819824371e-07, 2.27089458246136e-07, 2.65109018729386e-07, 3.09533787294144e-07,
>> 186 3.61450678951755e-07, 4.22132605720001e-07, 4.93070620012029e-07, 5.76011677884208e-07, 6.73003018378312e-07,
>> 187 7.86444334741884e-07, 9.19149125868381e-07, 1.07441686808177e-06, 1.25611794581916e-06, 1.46879363370693e-06,
>> 188 1.71777384258668e-06, 2.00931584092477e-06, 2.35076775595475e-06, 2.75076136411589e-06, 3.21943951980544e-06,
>> 189 3.76872457154335e-06, 4.41263530713955e-06, 5.16766139268237e-06, 6.05320597042372e-06, 7.09210911391032e-06,
>> 190 8.31126727283062e-06, 9.74236675732094e-06, 1.14227528119565e-05, 1.33964600351938e-05, 1.57154349593153e-05,
>> 191 1.84409877007598e-05, 2.16455169438847e-05, 2.54145614063097e-05, 2.98492416879238e-05, 3.50691694441802e-05,
>> 192 4.12159166620983e-05, 4.84571570931433e-05, 5.69916154060348e-05, 6.70549883575744e-05, 7.89270374851336e-05,
>> 193 9.29400960653008e-05, 1.09489286334512e-04, 1.29044808732337e-04, 1.52166746391861e-04, 1.79522929336019e-04,
>> 194 2.11910529073029e-04, 2.50282212267455e-04, 2.95777880652745e-04, 3.49763274771614e-04, 4.13877036477045e-04,
>> 195 4.90088229201730e-04, 5.80766832133355e-04, 6.88770389861787e-04, 8.17550860329037e-04, 9.71286825640136e-04,
>> 196 1.15504770108335e-03, 1.37499852012182e-03, 1.63865645831619e-03, 1.95521372859357e-03, 2.33594617839223e-03,
>> 197 2.79473334292108e-03, 3.34872458412049e-03, 4.01919834698563e-03, 4.83267910882034e-03, 5.82240174726978e-03,
>> 198 7.03024963447929e-03, 8.50934682352796e-03, 1.03275659802088e-02, 1.25723383030196e-02, 1.53573467148336e-02,
>> 199 1.88319961908362e-02, 2.31950693494851e-02, 2.87148467953825e-02, 3.57594981860148e-02, 4.48443278665925e-02,
>> 200 5.67077407341705e-02, 7.24383639141981e-02, 9.36982315562396e-02, 1.23138322807982e-01, 1.65231234966025e-01,
>> 201 2.28105620915395e-01, 3.28136673526381e-01, 5.03725917697149e-01, 8.70438998827130e-01, 2.00702876580162e+00
>> 202 };
>> 203
>> 204 const G4double G4GoudsmitSaundersonTable::fgSrcAValues[] = {
>> 205 1.04015860967600e+00, 1.04040150744807e+00, 1.04064741953810e+00, 1.04089640311900e+00, 1.04114851683194e+00,
>> 206 1.04140382083415e+00, 1.04166237684853e+00, 1.04192424821535e+00, 1.04218949994586e+00, 1.04245819877822e+00,
>> 207 1.04273041323563e+00, 1.04300621368681e+00, 1.04328567240902e+00, 1.04356886365371e+00, 1.04385586371482e+00,
>> 208 1.04414675100008e+00, 1.04444160610522e+00, 1.04474051189143e+00, 1.04504355356608e+00, 1.04535081876704e+00,
>> 209 1.04566239765056e+00, 1.04597838298311e+00, 1.04629887023725e+00, 1.04662395769178e+00, 1.04695374653644e+00,
>> 210 1.04728834098127e+00, 1.04762784837111e+00, 1.04797237930520e+00, 1.04832204776253e+00, 1.04867697123289e+00,
>> 211 1.04903727085420e+00, 1.04940307155639e+00, 1.04977450221209e+00, 1.05015169579471e+00, 1.05053478954418e+00,
>> 212 1.05092392514088e+00, 1.05131924888815e+00, 1.05172091190400e+00, 1.05212907032243e+00, 1.05254388550507e+00,
>> 213 1.05296552426367e+00, 1.05339415909403e+00, 1.05382996842239e+00, 1.05427313686456e+00, 1.05472385549904e+00,
>> 214 1.05518232215471e+00, 1.05564874171422e+00, 1.05612332643389e+00, 1.05660629628136e+00, 1.05709787929208e+00,
>> 215 1.05759831194585e+00, 1.05810783956480e+00, 1.05862671673423e+00, 1.05915520774789e+00, 1.05969358707925e+00,
>> 216 1.06024213988081e+00, 1.06080116251315e+00, 1.06137096310598e+00, 1.06195186215339e+00, 1.06254419314585e+00,
>> 217 1.06314830324155e+00, 1.06376455398007e+00, 1.06439332204139e+00, 1.06503500005380e+00, 1.06568999745437e+00,
>> 218 1.06635874140595e+00, 1.06704167777519e+00, 1.06773927217628e+00, 1.06845201108575e+00, 1.06918040303385e+00,
>> 219 1.06992497987890e+00, 1.07068629817131e+00, 1.07146494061474e+00, 1.07226151763259e+00, 1.07307666904867e+00,
>> 220 1.07391106589198e+00, 1.07476541233634e+00, 1.07564044778664e+00, 1.07653694912487e+00, 1.07745573313039e+00,
>> 221 1.07839765909007e+00, 1.07936363161610e+00, 1.08035460369071e+00, 1.08137157995935e+00, 1.08241562029607e+00,
>> 222 1.08348784366764e+00, 1.08458943232554e+00, 1.08572163635876e+00, 1.08688577864359e+00, 1.08808326023102e+00,
>> 223 1.08931556621724e+00, 1.09058427214773e+00, 1.09189105101191e+00, 1.09323768089211e+00, 1.09462605333836e+00,
>> 224 1.09605818254966e+00, 1.09753621545270e+00, 1.09906244278036e+00, 1.10063931126624e+00, 1.10226943708649e+00,
>> 225 1.10395562069822e+00, 1.10570086324418e+00, 1.10750838471683e+00, 1.10938164410283e+00, 1.11132436176008e+00,
>> 226 1.11334054431726e+00, 1.11543451242832e+00, 1.11761093176532e+00, 1.11987484769228e+00, 1.12223172413234e+00,
>> 227 1.12468748722310e+00, 1.12724857445230e+00, 1.12992199008194e+00, 1.13271536780716e+00, 1.13563704176072e+00,
>> 228 1.13869612717300e+00, 1.14190261223537e+00, 1.14526746300462e+00, 1.14880274353680e+00, 1.15252175386784e+00,
>> 229 1.15643918898390e+00, 1.16057132257199e+00, 1.16493622014373e+00, 1.16955398712386e+00, 1.17444705874537e+00,
>> 230 1.17964054016868e+00, 1.18516260723776e+00, 1.19104498083299e+00, 1.19732349105124e+00, 1.20403875167602e+00,
>> 231 1.21123697091984e+00, 1.21897093167774e+00, 1.22730118415524e+00, 1.23629750662007e+00, 1.24604070744941e+00,
>> 232 1.25662486545524e+00, 1.26816013838519e+00, 1.28077631555756e+00, 1.29462735591193e+00, 1.30989724673401e+00,
>> 233 1.32680765564328e+00, 1.34562805256148e+00, 1.36668928751138e+00, 1.39040208794754e+00, 1.41728269492652e+00,
>> 234 1.44798908275733e+00, 1.48337325073139e+00, 1.52455859525395e+00, 1.57305765486728e+00, 1.63095721573069e+00,
>> 235 1.70122060331290e+00, 1.78820418157127e+00, 1.89858943787341e+00, 2.04318263724065e+00, 2.24070141295433e+00,
>> 236 2.52670054341548e+00, 2.97823240973856e+00, 3.80070470423559e+00, 5.80504410098720e+00, 0.0
>> 237 };
>> 238
>> 239 const G4double G4GoudsmitSaundersonTable::fgSrcBValues[] = {
>> 240 4.85267101555493e-02, 4.87971711674324e-02, 4.90707875373474e-02, 4.93476163203902e-02, 4.96277159833787e-02,
>> 241 4.99111464494480e-02, 5.01979691443985e-02, 5.04882470448502e-02, 5.07820447282697e-02, 5.10794284249799e-02,
>> 242 5.13804660722384e-02, 5.16852273704654e-02, 5.19937838417706e-02, 5.23062088908123e-02, 5.26225778681911e-02,
>> 243 5.29429681364105e-02, 5.32674591386191e-02, 5.35961324701870e-02, 5.39290719533364e-02, 5.42663637149114e-02,
>> 244 5.46080962674657e-02, 5.49543605938995e-02, 5.53052502357015e-02, 5.56608613851103e-02, 5.60212929813128e-02,
>> 245 5.63866468109177e-02, 5.67570276129776e-02, 5.71325431886598e-02, 5.75133045160478e-02, 5.78994258700939e-02,
>> 246 5.82910249481984e-02, 5.86882230016298e-02, 5.90911449731285e-02, 5.94999196410672e-02, 5.99146797704773e-02,
>> 247 6.03355622714283e-02, 6.07627083650703e-02, 6.11962637578703e-02, 6.16363788244832e-02, 6.20832087997576e-02,
>> 248 6.25369139805019e-02, 6.29976599373900e-02, 6.34656177379503e-02, 6.39409641809610e-02, 6.44238820432201e-02,
>> 249 6.49145603393270e-02, 6.54131945953447e-02, 6.59199871371941e-02, 6.64351473947443e-02, 6.69588922226319e-02,
>> 250 6.74914462388500e-02, 6.80330421823362e-02, 6.85839212907651e-02, 6.91443337000157e-02, 6.97145388666155e-02,
>> 251 7.02948060148960e-02, 7.08854146105257e-02, 7.14866548622191e-02, 7.20988282536816e-02, 7.27222481079186e-02,
>> 252 7.33572401862953e-02, 7.40041433248063e-02, 7.46633101103861e-02, 7.53351076002229e-02, 7.60199180872846e-02,
>> 253 7.67181399156741e-02, 7.74301883495570e-02, 7.81564964999089e-02, 7.88975163136540e-02, 7.96537196301176e-02,
>> 254 8.04255993102895e-02, 8.12136704447979e-02, 8.20184716471344e-02, 8.28405664392583e-02, 8.36805447373290e-02,
>> 255 8.45390244462388e-02, 8.54166531722937e-02, 8.63141100644354e-02, 8.72321077953886e-02, 8.81713946953474e-02,
>> 256 8.91327570520271e-02, 9.01170215924217e-02, 9.11250581632225e-02, 9.21577826286684e-02, 9.32161600065677e-02,
>> 257 9.43012078656787e-02, 9.54140000100007e-02, 9.65556704785876e-02, 9.77274178926762e-02, 9.89305101855936e-02,
>> 258 1.00166289755146e-01, 1.01436179082793e-01, 1.02741686869338e-01, 1.04084414742952e-01, 1.05466064602181e-01,
>> 259 1.06888446664513e-01, 1.08353488300132e-01, 1.09863243740694e-01, 1.11419904764858e-01, 1.13025812475804e-01,
>> 260 1.14683470301675e-01, 1.16395558367910e-01, 1.18164949411119e-01, 1.19994726428692e-01, 1.21888202285954e-01,
>> 261 1.23848941535858e-01, 1.25880784744097e-01, 1.27987875657392e-01, 1.30174691605324e-01, 1.32446077587840e-01,
>> 262 1.34807284573851e-01, 1.37264012622799e-01, 1.39822459544273e-01, 1.42489375933601e-01, 1.45272127568431e-01,
>> 263 1.48178766328292e-01, 1.51218111012368e-01, 1.54399839689173e-01, 1.57734595526106e-01, 1.61234108430976e-01,
>> 264 1.64911335308903e-01, 1.68780622319388e-01, 1.72857893239124e-01, 1.77160868934300e-01, 1.81709324071779e-01,
>> 265 1.86525388617776e-01, 1.91633903472561e-01, 1.97062841888240e-01, 2.02843811271485e-01, 2.09012653799732e-01,
>> 266 2.15610169273952e-01, 2.22682990203222e-01, 2.30284647840434e-01, 2.38476879578463e-01, 2.47331243935506e-01,
>> 267 2.56931130997193e-01, 2.67374286122045e-01, 2.78776006654390e-01, 2.91273230921782e-01, 3.05029824532737e-01,
>> 268 3.20243494424876e-01, 3.37154947792243e-01, 3.56060196110919e-01, 3.77327342746089e-01, 4.01419886817026e-01,
>> 269 4.28929703940887e-01, 4.60624750236354e-01, 4.97519791888632e-01, 5.40984294142651e-01, 5.92912498016053e-01,
>> 270 6.56002088728957e-01, 7.34232298458703e-01, 8.33731313414106e-01, 9.64463147693015e-01, 1.14381339640412e+00,
>> 271 1.40517335796248e+00, 1.82227062151045e+00, 2.59912058457152e+00, 4.62773895951686e+00, 0.0
>> 272 };
>> 273
>> 274
>> 275 G4double G4GoudsmitSaundersonTable::fgInverseQ2CDFs[fgNumLambdas*fgNumLamG1*fgNumUvalues] = {0.0};
>> 276 G4double G4GoudsmitSaundersonTable::fgInterParamsA2[fgNumLambdas*fgNumLamG1*fgNumUvalues] = {0.0};
>> 277 G4double G4GoudsmitSaundersonTable::fgInterParamsB2[fgNumLambdas*fgNumLamG1*fgNumUvalues] = {0.0};
>> 278 G4double G4GoudsmitSaundersonTable::fgInverseQ2CDFsII[fgNumLambdas*fgNumLamG1II*fgNumUvalues] = {0.0};
>> 279 G4double G4GoudsmitSaundersonTable::fgInterParamsA2II[fgNumLambdas*fgNumLamG1II*fgNumUvalues] = {0.0};
>> 280 G4double G4GoudsmitSaundersonTable::fgInterParamsB2II[fgNumLambdas*fgNumLamG1II*fgNumUvalues] = {0.0};
>> 281
>> 282 std::vector<G4double>* G4GoudsmitSaundersonTable::fgMoliereBc = 0;
>> 283 std::vector<G4double>* G4GoudsmitSaundersonTable::fgMoliereXc2 = 0;
>> 284
>> 285 G4bool G4GoudsmitSaundersonTable::fgIsInitialised = FALSE;
>> 286
>> 287 ////////////////////////////////////////////////////////////////////////////////
>> 288 void G4GoudsmitSaundersonTable::Initialise(){
>> 289 if(!fgIsInitialised){
>> 290 LoadMSCData();
>> 291 LoadMSCDataII();
>> 292 fgIsInitialised = TRUE;
>> 293 }
111 294
112 G4GoudsmitSaundersonTable::G4GoudsmitSaunderso << 295 InitMoliereMSCParams();
113 fIsElectron = iselectron; <<
114 // set initial values: final values will be <<
115 fLogLambda0 = 0.; // wi <<
116 fLogDeltaLambda = 0.; // wi <<
117 fInvLogDeltaLambda = 0.; // wi <<
118 fInvDeltaQ1 = 0.; // wi <<
119 fDeltaQ2 = 0.; // wi <<
120 fInvDeltaQ2 = 0.; // wi <<
121 // <<
122 fLowEnergyLimit = 0.1*CLHEP::keV; // w <<
123 fHighEnergyLimit = 100.0*CLHEP::MeV; // w <<
124 // <<
125 fIsMottCorrection = false; // w <<
126 fIsPWACorrection = false; // w <<
127 fMottCorrection = nullptr; <<
128 // <<
129 fNumSPCEbinPerDec = 3; <<
130 } 296 }
131 297
132 G4GoudsmitSaundersonTable::~G4GoudsmitSaunders << 298 ////////////////////////////////////////////////////////////////////////////////
133 for (std::size_t i=0; i<gGSMSCAngularDistrib << 299 G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable(){
134 if (gGSMSCAngularDistributions1[i]) { << 300 if(fgMoliereBc){
135 delete [] gGSMSCAngularDistributions1[i] << 301 delete fgMoliereBc;
136 delete [] gGSMSCAngularDistributions1[i] << 302 fgMoliereBc = 0;
137 delete [] gGSMSCAngularDistributions1[i] << 303 }
138 delete gGSMSCAngularDistributions1[i]; << 304 if(fgMoliereXc2){
139 } << 305 delete fgMoliereXc2;
140 } << 306 fgMoliereXc2 = 0;
141 gGSMSCAngularDistributions1.clear(); << 307 }
142 for (std::size_t i=0; i<gGSMSCAngularDistrib << 308 fgIsInitialised = FALSE;
143 if (gGSMSCAngularDistributions2[i]) { << 309 }
144 delete [] gGSMSCAngularDistributions2[i] << 310
145 delete [] gGSMSCAngularDistributions2[i] << 311 ////////////////////////////////////////////////////////////////////////////////
146 delete [] gGSMSCAngularDistributions2[i] << 312 // sample cos(theta) i.e. angular deflection from the precomputed angular
147 delete gGSMSCAngularDistributions2[i]; << 313 // distributions in the real multiple scattering case
148 } << 314 G4double G4GoudsmitSaundersonTable::SampleCosTheta(G4double lambdavalue, G4double lamG1value,
149 } << 315 G4double screeningparam, G4double rndm1, G4double rndm2,
150 gGSMSCAngularDistributions2.clear(); << 316 G4double rndm3){
151 if (fMottCorrection) { << 317 const G4double lnl0 = 0.0; //log(fgLambdaValues[0]);
152 delete fMottCorrection; << 318 const G4double invlnl = 6.5144172285e+00; //1./log(fgLambdaValues[i+1]/fgLambdaValues[i])
153 fMottCorrection = nullptr; << 319 G4double logLambdaVal = G4Log(lambdavalue);
154 } << 320 G4int lambdaindx = (G4int)((logLambdaVal-lnl0)*invlnl);
155 // clear scp correction data << 321
156 for (std::size_t imc=0; imc<fSCPCPerMatCuts. << 322 const G4double dd = 1.0/0.0499; // delta
157 if (fSCPCPerMatCuts[imc]) { << 323 G4int lamg1indx = (G4int)((lamG1value-fgLamG1Values[0])*dd);
158 fSCPCPerMatCuts[imc]->fVSCPC.clear(); << 324
159 delete fSCPCPerMatCuts[imc]; << 325 G4double probOfLamJ = G4Log(fgLambdaValues[lambdaindx+1]/lambdavalue)*invlnl;
160 } << 326 if(rndm1 > probOfLamJ)
161 } << 327 ++lambdaindx;
162 fSCPCPerMatCuts.clear(); << 328
163 // << 329 // store 1.0/(fgLamG1Values[lamg1indx+1]-fgLamG1Values[lamg1indx]) const value
164 gIsInitialised = false; << 330 const G4double onePerLamG1Diff = dd;
165 } << 331 G4double probOfLamG1J = (fgLamG1Values[lamg1indx+1]-lamG1value)*onePerLamG1Diff;
>> 332 if(rndm2 > probOfLamG1J)
>> 333 ++lamg1indx;
>> 334
>> 335 G4int begin = lambdaindx*(fgNumLamG1*fgNumUvalues)+lamg1indx*fgNumUvalues;
>> 336
>> 337 // sample bin
>> 338 G4int indx = (G4int)(rndm3*100); // value/delatU ahol deltaU = 0.01; find u-indx
>> 339 G4double cp1 = fgUValues[indx];
>> 340 // G4double cp2 = fgUValues[indx+1];
>> 341 indx +=begin;
>> 342
>> 343 G4double u1 = fgInverseQ2CDFs[indx];
>> 344 G4double u2 = fgInverseQ2CDFs[indx+1];
>> 345
>> 346 G4double dum1 = rndm3 - cp1;
>> 347 const G4double dum2 = 0.01;//cp2-cp1; // this is constans 0.01
>> 348 const G4double dum22 = 0.0001;
>> 349 // rational interpolation of the inverse CDF
>> 350 G4double sample= u1+(1.0+fgInterParamsA2[indx]+fgInterParamsB2[indx])*dum1*dum2/(dum22+
>> 351 dum1*dum2*fgInterParamsA2[indx]+fgInterParamsB2[indx]*dum1*dum1)*(u2-u1);
166 352
167 void G4GoudsmitSaundersonTable::Initialise(G4d << 353 // transform back u to cos(theta) :
168 fLowEnergyLimit = lownergylimit; << 354 // -compute parameter value a = w2*screening_parameter
169 fHighEnergyLimit = highenergylimit; << 355 G4double t = logLambdaVal;
170 G4double lLambdaMin = G4Log(gLAMBMIN); << 356 G4double dum;
171 G4double lLambdaMax = G4Log(gLAMBMAX); << 357 if(lambdavalue >= 10.0)
172 fLogLambda0 = lLambdaMin; << 358 dum = 0.5*(-2.77164+t*( 2.94874-t*(0.1535754-t*0.00552888) ));
173 fLogDeltaLambda = (lLambdaMax-lLambdaMin << 359 else
174 fInvLogDeltaLambda = 1./fLogDeltaLambda; << 360 dum = 0.5*(1.347+t*(0.209364-t*(0.45525-t*(0.50142-t*0.081234))));
175 fInvDeltaQ1 = 1./((gQMAX1-gQMIN1)/(g << 361
176 fDeltaQ2 = (gQMAX2-gQMIN2)/(gQNUM << 362 G4double a = dum*(lambdavalue+4.0)*screeningparam;
177 fInvDeltaQ2 = 1./fDeltaQ2; << 363
178 // load precomputed angular distributions an << 364 // this is the sampled cos(theta)
179 // these are particle independet => they go << 365 return (2.0*a*sample)/(1.0-sample+a);
180 if (!gIsInitialised) { << 366 }
181 // load pre-computed GS angular distributi << 367
182 LoadMSCData(); << 368 ////////////////////////////////////////////////////////////////////////////////
183 gIsInitialised = true; << 369 // sample cos(theta) i.e. angular deflection from the precomputed angular
184 } << 370 // distributions in the real multiple scattering case (For the second grid)
185 InitMoliereMSCParams(); << 371 G4double G4GoudsmitSaundersonTable::SampleCosThetaII(G4double lambdavalue, G4double lamG1value,
186 // Mott-correction: particle(e- or e+) depen << 372 G4double screeningparam, G4double rndm1, G4double rndm2,
187 if (fIsMottCorrection) { << 373 G4double rndm3){
188 if (!fMottCorrection) { << 374 const G4double lnl0 = 0.0; //log(fgLambdaValues[0]);
189 fMottCorrection = new G4GSMottCorrection << 375 const G4double invlnl = 6.5144172285e+00; //1./log(fgLambdaValues[i+1]/fgLambdaValues[i])
190 } << 376 G4double logLambdaVal = G4Log(lambdavalue);
191 fMottCorrection->Initialise(); << 377 G4int lambdaindx = (G4int)((logLambdaVal-lnl0)*invlnl);
192 } << 378
193 // init scattering power correction data; us << 379 const G4double dd = 1.0/0.333; // delta
194 // (Moliere's parameters must be initialised << 380 G4int lamg1indx = (G4int)((lamG1value-fgLamG1ValuesII[0])*dd);
195 if (fMottCorrection) { << 381
196 InitSCPCorrection(); << 382 G4double probOfLamJ = G4Log(fgLambdaValues[lambdaindx+1]/lambdavalue)*invlnl;
197 } << 383 if(rndm1 > probOfLamJ)
198 } << 384 ++lambdaindx;
>> 385
>> 386 // store 1.0/(fgLamG1Values[lamg1indx+1]-fgLamG1Values[lamg1indx]) const value
>> 387 const G4double onePerLamG1Diff = dd;
>> 388 G4double probOfLamG1J = (fgLamG1ValuesII[lamg1indx+1]-lamG1value)*onePerLamG1Diff;
>> 389 if(rndm2 > probOfLamG1J)
>> 390 ++lamg1indx;
>> 391
>> 392 // protection against cases when the sampled lamG1 values are not possible due
>> 393 // to the fact that G1<=1 (G1 = lam_el/lam_tr1)
>> 394 // -- in theory it should never happen when lamg1indx=0 but check it just to be sure
>> 395 if(fgLamG1ValuesII[lamg1indx]>lambdavalue && lamg1indx>0)
>> 396 --lamg1indx;
>> 397
>> 398
>> 399 G4int begin = lambdaindx*(fgNumLamG1II*fgNumUvalues)+lamg1indx*fgNumUvalues;
>> 400
>> 401 // sample bin
>> 402 G4int indx = (G4int)(rndm3*100); // value/delatU ahol deltaU = 0.01; find u-indx
>> 403 G4double cp1 = fgUValues[indx];
>> 404 // G4double cp2 = fgUValues[indx+1];
>> 405 indx +=begin;
>> 406
>> 407 G4double u1 = fgInverseQ2CDFsII[indx];
>> 408 G4double u2 = fgInverseQ2CDFsII[indx+1];
>> 409
>> 410 G4double dum1 = rndm3 - cp1;
>> 411 const G4double dum2 = 0.01;//cp2-cp1; // this is constans 0.01
>> 412 const G4double dum22 = 0.0001;
>> 413 // rational interpolation of the inverse CDF
>> 414 G4double sample= u1+(1.0+fgInterParamsA2II[indx]+fgInterParamsB2II[indx])*dum1*dum2/(dum22+
>> 415 dum1*dum2*fgInterParamsA2II[indx]+fgInterParamsB2II[indx]*dum1*dum1)*(u2-u1);
>> 416
>> 417 // transform back u to cos(theta) :
>> 418 // -compute parameter value a = w2*screening_parameter
>> 419 G4double t = logLambdaVal;
>> 420 G4double dum = 0.;
>> 421 if(lambdavalue >= 10.0)
>> 422 dum = 0.5*(-2.77164+t*( 2.94874-t*(0.1535754-t*0.00552888) ));
>> 423 else
>> 424 dum = 0.5*(1.347+t*(0.209364-t*(0.45525-t*(0.50142-t*0.081234))));
199 425
>> 426 G4double a = dum*(lambdavalue+4.0)*screeningparam;
200 427
>> 428 // this is the sampled cos(theta)
>> 429 return (2.0*a*sample)/(1.0-sample+a);
>> 430 }
>> 431
>> 432 ////////////////////////////////////////////////////////////////////////////////
201 // samplig multiple scattering angles cos(thet 433 // samplig multiple scattering angles cos(theta) and sin(thata)
202 // - including no-scattering, single, "few" s << 434 // including no-scattering, single, "few" scattering cases as well
203 // - Mott-correction will be included if it w << 435 void G4GoudsmitSaundersonTable::Sampling(G4double lambdavalue, G4double lamG1value,
204 // lambdaval : s/lambda_el << 436 G4double scrPar, G4double &cost, G4double &sint){
205 // qval : s/lambda_el G1 <<
206 // scra : screening parameter <<
207 // cost : will be the smapled cos(theta) <<
208 // sint : will be the smapled sin(theta) <<
209 // lekin : logarithm of the current kineti <<
210 // beta2 : the corresponding beta square <<
211 // matindx : index of the current material <<
212 // returns true if it was msc <<
213 G4bool G4GoudsmitSaundersonTable::Sampling(G4d <<
214 G4d <<
215 GSM <<
216 G4d <<
217 G4double rand0 = G4UniformRand(); 437 G4double rand0 = G4UniformRand();
218 G4double expn = G4Exp(-lambdaval); << 438 G4double expn = G4Exp(-lambdavalue);
>> 439
219 // 440 //
220 // no scattering case 441 // no scattering case
221 if (rand0<expn) { << 442 //
>> 443 if(rand0 < expn) {
222 cost = 1.0; 444 cost = 1.0;
223 sint = 0.0; 445 sint = 0.0;
224 return false; << 446 return;
225 } 447 }
>> 448
>> 449 //
>> 450 // single scattering case : sample (analitically) from the single scattering PDF
>> 451 // that corresponds to the modified-Wentzel DCS i.e.
>> 452 // Wentzel DCS that gives back the PWA first-transport mfp
226 // 453 //
227 // single scattering case : sample from the << 454 if( rand0 < (1.+lambdavalue)*expn) {
228 // - Mott-correction will be included if it << 455 G4double rand1 = G4UniformRand();
229 if (rand0<(1.+lambdaval)*expn) { << 456 // sampling 1-cos(theta)
230 // cost is sampled in SingleScattering() << 457 G4double dum0 = 2.0*scrPar*rand1/(1.0 - rand1 + scrPar);
231 cost = SingleScattering(lambdaval, scra, l <<
232 // add protections 458 // add protections
233 if (cost<-1.0) cost = -1.0; << 459 if(dum0 < 0.0) dum0 = 0.0;
234 if (cost>1.0) cost = 1.0; << 460 else if(dum0 > 2.0 ) dum0 = 2.0;
235 // compute sin(theta) from the sampled cos << 461 // compute cos(theta) and sin(theta)
236 G4double dum0 = 1.-cost; << 462 cost = 1.0 - dum0;
237 sint = std::sqrt(dum0*(2.0-dum0)); << 463 sint = std::sqrt(dum0*(2.0 - dum0));
238 return false; << 464 return;
239 } 465 }
>> 466
240 // 467 //
241 // handle this case: 468 // handle this case:
242 // -lambdaval < 1 i.e. mean #elastic ev << 469 // -lambdavalue < 1 i.e. mean #elastic events along the step is < 1 but
243 // the currently sampled case is not 0 470 // the currently sampled case is not 0 or 1 scattering. [Our minimal
244 // lambdaval (that we have precomputed << 471 // lambdavalue (that we have precomputed, transformed angular distributions
245 // stored in a form of equally probabe 472 // stored in a form of equally probabe intervalls together with rational
246 // interp. parameters) is 1.] 473 // interp. parameters) is 1.]
247 // -probability of having n elastic eve 474 // -probability of having n elastic events follows Poisson stat. with
248 // lambdaval parameter. << 475 // lambdavalue parameter.
249 // -the max. probability (when lambdava << 476 // -the max. probability (when lambdavalue=1) of having more than one
250 // elastic events is 0.2642411 and the 477 // elastic events is 0.2642411 and the prob of having 2,3,..,n elastic
251 // events decays rapidly with n. So se 478 // events decays rapidly with n. So set a max n to 10.
252 // -sampling of this cases is done in a 479 // -sampling of this cases is done in a one-by-one single elastic event way
253 // where the current #elastic event is 480 // where the current #elastic event is sampled from the Poisson distr.
254 if (lambdaval<1.0) { << 481 // -(other possibility would be to use lambdavalue=1 distribution because
255 G4double prob, cumprob; << 482 // they are very close)
256 prob = cumprob = expn; << 483 if(lambdavalue < 1.0) {
257 G4double curcost,cursint; << 484 G4double prob, cumprob;
258 // init cos(theta) and sin(theta) to the z << 485 prob = cumprob = expn;
259 cost = 1.0; << 486 G4double curcost,cursint;
260 sint = 0.0; << 487 // init cos(theta) and sin(theta) to the zero scattering values
261 for (G4int iel=1; iel<10; ++iel) { << 488 cost = 1.0;
262 // prob of having iel scattering from Po << 489 sint = 0.0;
263 prob *= lambdaval/(G4double)iel; << 490
264 cumprob += prob; << 491 for(G4int iel = 1; iel < 10; ++iel){
265 // << 492 // prob of having iel scattering from Poisson
266 //sample cos(theta) from the singe scatt << 493 prob *= lambdavalue/(G4double)iel;
267 // - Mott-correction will be included if << 494 cumprob += prob;
268 curcost = SingleScattering(lambdav << 495 //
269 G4double dum0 = 1.-curcost; << 496 //sample cos(theta) from the singe scattering pdf
270 cursint = dum0*(2.0-dum0); // sin^ << 497 //
271 // << 498 G4double rand1 = G4UniformRand();
272 // if we got current deflection that is << 499 // sampling 1-cos(theta)
273 // then update cos(theta) sin(theta) << 500 G4double dum0 = 2.0*scrPar*rand1/(1.0 - rand1 + scrPar);
274 if (cursint>1.0e-20) { << 501 // compute cos(theta) and sin(theta)
275 cursint = std::sqrt(cursint); << 502 curcost = 1.0 - dum0;
276 G4double curphi = CLHEP::twopi*G4Unifo << 503 cursint = dum0*(2.0 - dum0);
277 cost = cost*curcost-sint*cu << 504 //
278 sint = std::sqrt(std::max(0 << 505 // if we got current deflection that is not too small
279 } << 506 // then update cos(theta) sin(theta)
280 // << 507 if(cursint > 1.0e-20){
281 // check if we have done enough scatteri << 508 cursint = std::sqrt(cursint);
282 if (rand0<cumprob) { << 509 G4double curphi = CLHEP::twopi*G4UniformRand();
283 return false; << 510 cost = cost*curcost - sint*cursint*std::cos(curphi);
284 } << 511 sint = std::sqrt(std::max(0.0, (1.0-cost)*(1.0+cost)));
285 } << 512 }
286 // if reached the max iter i.e. 10 << 513 //
287 return false; << 514 // check if we have done enough scattering i.e. sampling from the Poisson
>> 515 if( rand0 < cumprob)
>> 516 return;
>> 517 }
>> 518 // if reached the max iter i.e. 10
>> 519 return;
288 } 520 }
>> 521
>> 522
>> 523 /// **** let it izotropic if we are above the grid i.e. true path length is long
>> 524 // IT IS DONE NOW IN THE MODEL
>> 525 //if(lamG1value>fgLamG1Values[fgNumLamG1-1]) {
>> 526 // cost = 1.-2.*G4UniformRand();
>> 527 // sint = std::sqrt((1.-cost)*(1.+cost));
>> 528 // return;
>> 529 //}
>> 530
>> 531
289 // 532 //
290 // multiple scattering case with lambdavalue 533 // multiple scattering case with lambdavalue >= 1:
291 // - use the precomputed and transformed G 534 // - use the precomputed and transformed Goudsmit-Saunderson angular
292 // distributions to sample cos(theta) << 535 // distributions that are stored in a form of equally probabe intervalls
293 // - Mott-correction will be included if i << 536 // together with the corresponding rational interp. parameters
294 cost = SampleCosTheta(lambdaval, qval, scra, << 537 //
295 // add protections 538 // add protections
296 if (cost<-1.0) cost = -1.0; << 539 if(lamG1value<fgLamG1Values[0]) lamG1value = fgLamG1Values[0];
297 if (cost> 1.0) cost = 1.0; << 540 if(lamG1value>fgLamG1ValuesII[fgNumLamG1II-1]) lamG1value = fgLamG1ValuesII[fgNumLamG1II-1];
298 // compute cos(theta) and sin(theta) from th <<
299 G4double dum0 = 1.0-cost; <<
300 sint = std::sqrt(dum0*(2.0-dum0)); <<
301 // return true if it was msc <<
302 return true; <<
303 } <<
304 541
>> 542 // sample 1-cos(theta) :: depending on the grids
>> 543 G4double dum0 = 0.0;
>> 544 if(lamG1value<fgLamG1Values[fgNumLamG1-1]) {
>> 545 dum0 = SampleCosTheta(lambdavalue, lamG1value, scrPar, G4UniformRand(),
>> 546 G4UniformRand(), G4UniformRand());
305 547
306 G4double G4GoudsmitSaundersonTable::SampleCosT << 548 } else {
307 << 549 dum0 = SampleCosThetaII(lambdavalue, lamG1value, scrPar, G4UniformRand(),
308 << 550 G4UniformRand(), G4UniformRand());
309 <<
310 G4double cost = 1.; <<
311 // determine the base GS angular distributio <<
312 if (isfirst) { <<
313 *gsDtr = GetGSAngularDtr(scra, lambdaval, <<
314 } <<
315 // sample cost from the GS angular distribut <<
316 cost = SampleGSSRCosTheta(*gsDtr, transfPar) <<
317 // Mott-correction if it was requested by th <<
318 if (fIsMottCorrection && *gsDtr) { // no <<
319 static const G4int nlooplim = 1000; <<
320 G4int nloop = 0 ; // rejection loop <<
321 // G4int ekindx = -1; // evaluate only <<
322 // G4int deltindx = -1 ; // evaluate onl <<
323 G4double val = fMottCorrection->GetMo <<
324 while (G4UniformRand()>val && ++nloop<nloo <<
325 // sampling cos(theta) <<
326 cost = SampleGSSRCosTheta(*gsDtr, transf <<
327 val = fMottCorrection->GetMottRejection <<
328 }; <<
329 } 551 }
330 return cost; <<
331 } <<
332 552
>> 553 // add protections
>> 554 if(dum0 < 0.0) dum0 = 0.0;
>> 555 if(dum0 > 2.0 ) dum0 = 2.0;
333 556
334 // returns with cost sampled from the GS angul << 557 // compute cos(theta) and sin(theta)
335 G4double G4GoudsmitSaundersonTable::SampleGSSR << 558 cost = 1.0-dum0;
336 // check if isotropic theta (i.e. cost is un << 559 sint = std::sqrt(dum0*(2.0 - dum0));
337 if (!gsDtr) { << 560 }
338 return 1.-2.0*G4UniformRand(); << 561
339 } << 562 ////////////////////////////////////////////////////////////////////////////////
340 // << 563 G4double G4GoudsmitSaundersonTable::GetScreeningParam(G4double G1){
341 // sampling form the selected distribution << 564 if(G1<fgG1Values[0]) return fgScreeningParam[0];
342 G4double ndatm1 = gsDtr->fNumData-1.; << 565 if(G1>fgG1Values[fgNumScreeningParams-1]) return fgScreeningParam[fgNumScreeningParams-1];
343 G4double delta = 1.0/ndatm1; << 566 // normal case
344 // determine lower cumulative bin inidex << 567 const G4double lng10 = -2.30258509299405e+01; //log(fgG1Values[0]);
345 G4double rndm = G4UniformRand(); << 568 const G4double invlng1 = 6.948711710452e+00; //1./log(fgG1Values[i+1]/fgG1Values[i])
346 G4int indxl = rndm*ndatm1; << 569 G4double logG1 = G4Log(G1);
347 G4double aval = rndm-indxl*delta; << 570 G4int k = (G4int)((logG1-lng10)*invlng1);
348 G4double dum0 = delta*aval; << 571 return G4Exp(logG1*fgSrcAValues[k])*fgSrcBValues[k]; // log-log linear
349 << 572 }
350 G4double dum1 = (1.0+gsDtr->fParamA[indxl] << 573
351 G4double dum2 = delta*delta + gsDtr->fPara << 574 ////////////////////////////////////////////////////////////////////////////////
352 G4double sample = gsDtr->fUValues[indxl] + << 575 void G4GoudsmitSaundersonTable::LoadMSCData(){
353 // transform back u to cos(theta) : << 576 G4double dum;
354 // this is the sampled cos(theta) = (2.0*par << 577 char basename[512];
355 return 1.-(2.0*transfpar*sample)/(1.0-sample << 578 char* path = getenv("G4LEDATA");
356 } << 579 if (!path) {
>> 580 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
>> 581 FatalException,
>> 582 "Environment variable G4LEDATA not defined");
>> 583 return;
>> 584 }
>> 585
>> 586 std::string pathString(path);
>> 587 sprintf(basename,"inverseq2CDF");
>> 588 for(int k = 0; k < fgNumLambdas; ++k){
>> 589 char fname[512];
>> 590 sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
>> 591 std::ifstream infile(fname,std::ios::in);
>> 592 if(!infile.is_open()){
>> 593 char msgc[512];
>> 594 sprintf(msgc,"Cannot open file: %s .",fname);
>> 595 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
>> 596 FatalException,
>> 597 msgc);
>> 598 return;
>> 599 }
357 600
>> 601 int limk = k*fgNumUvalues*fgNumLamG1;
>> 602 for(int i = 0; i < fgNumUvalues; ++i)
>> 603 for(int j = 0; j < fgNumLamG1+1; ++j)
>> 604 if(j==0) infile >> dum;
>> 605 else infile >> fgInverseQ2CDFs[limk+(j-1)*fgNumUvalues+i];
>> 606
>> 607 infile.close();
>> 608 }
>> 609
>> 610
>> 611 sprintf(basename,"inverseq2CDFRatInterpPA");
>> 612 for(int k = 0; k < fgNumLambdas; ++k){
>> 613 char fname[512];
>> 614 sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
>> 615 std::ifstream infile(fname,std::ios::in);
>> 616 if(!infile.is_open()){
>> 617 char msgc[512];
>> 618 sprintf(msgc,"Cannot open file: %s .",fname);
>> 619 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
>> 620 FatalException,
>> 621 msgc);
>> 622 return;
>> 623 }
358 624
359 // determine the GS angular distribution we ne << 625 int limk = k*fgNumUvalues*fgNumLamG1;
360 G4GoudsmitSaundersonTable::GSMSCAngularDtr* G4 << 626 for(int i = 0; i < fgNumUvalues; ++i)
361 << 627 for(int j = 0; j < fgNumLamG1+1; ++j)
362 GSMSCAngularDtr *dtr = nullptr; << 628 if(j==0) infile >> dum;
363 G4bool first = false; << 629 else infile >> fgInterParamsA2[limk+(j-1)*fgNumUvalues+i];
364 // isotropic cost above gQMAX2 (i.e. dtr sta << 630
365 if (qval<gQMAX2) { << 631 infile.close();
366 G4int lamIndx = -1; // lambda value in << 632 }
367 G4int qIndx = -1; // lambda value in << 633
368 // init to second grid Q values << 634 sprintf(basename,"inverseq2CDFRatInterpPB");
369 G4int numQVal = gQNUM2; << 635 for(int k = 0; k < fgNumLambdas; ++k){
370 G4double minQVal = gQMIN2; << 636 char fname[512];
371 G4double invDelQ = fInvDeltaQ2; << 637 sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
372 G4double pIndxH = 0.; // probability of << 638 std::ifstream infile(fname,std::ios::in);
373 // check if first or second grid needs to << 639 if(!infile.is_open()){
374 if (qval<gQMIN2) { // first grid << 640 char msgc[512];
375 first = true; << 641 sprintf(msgc,"Cannot open file: %s .",fname);
376 // protect against qval<gQMIN1 << 642 G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
377 if (qval<gQMIN1) { << 643 FatalException,
378 qval = gQMIN1; << 644 msgc);
379 qIndx = 0; << 645 return;
380 //pIndxH = 0.; << 646 }
381 } << 647
382 // set to first grid Q values << 648 int limk = k*fgNumUvalues*fgNumLamG1;
383 numQVal = gQNUM1; << 649 for(int i = 0; i < fgNumUvalues; ++i)
384 minQVal = gQMIN1; << 650 for(int j = 0; j < fgNumLamG1+1; ++j)
385 invDelQ = fInvDeltaQ1; << 651 if(j==0) infile >> dum;
386 } << 652 else infile >> fgInterParamsB2[limk+(j-1)*fgNumUvalues+i];
387 // make sure that lambda = s/lambda_el is << 653
388 // lambda<gLAMBMIN=1 is already handeled b << 654 infile.close();
389 if (lambdaval>=gLAMBMAX) { << 655 }
390 lambdaval = gLAMBMAX-1.e-8; << 656 }
391 lamIndx = gLAMBNUM-1; << 657
392 } << 658 ////////////////////////////////////////////////////////////////////////////////
393 G4double lLambda = G4Log(lambdaval); << 659 // Load the second grid data
394 // << 660 ////////////////////////////////////////////////////////////////////////////////
395 // determine lower lambda (=s/lambda_el) i << 661 void G4GoudsmitSaundersonTable::LoadMSCDataII(){
396 if (lamIndx<0) { << 662 G4double dum;
397 pIndxH = (lLambda-fLogLambda0)*fInvLogD << 663 char basename[512];
398 lamIndx = (G4int)(pIndxH); // low << 664 char* path = getenv("G4LEDATA");
399 pIndxH = pIndxH-lamIndx; // proba << 665 if (!path) {
400 if (G4UniformRand()<pIndxH) { << 666 G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
401 ++lamIndx; << 667 FatalException,
402 } << 668 "Environment variable G4LEDATA not defined");
403 } << 669 return;
404 // <<
405 // determine lower Q (=s/lambda_el G1) ind <<
406 if (qIndx<0) { <<
407 pIndxH = (qval-minQVal)*invDelQ; <<
408 qIndx = (G4int)(pIndxH); // low <<
409 pIndxH = pIndxH-qIndx; <<
410 if (G4UniformRand()<pIndxH) { <<
411 ++qIndx; <<
412 } <<
413 } <<
414 // set indx <<
415 G4int indx = lamIndx*numQVal+qIndx; <<
416 if (first) { <<
417 dtr = gGSMSCAngularDistributions1[indx]; <<
418 } else { <<
419 dtr = gGSMSCAngularDistributions2[indx]; <<
420 } <<
421 // dtr might be nullptr that indicates iso <<
422 // - if the selected lamIndx, qIndx corres <<
423 // G1 should always be < 1 and if G1 is <<
424 // <<
425 // compute the transformation parameter <<
426 if (lambdaval>10.0) { <<
427 transfpar = 0.5*(-2.77164+lLambda*( 2.94 <<
428 } else { <<
429 transfpar = 0.5*(1.347+lLambda*(0.209364 <<
430 } <<
431 transfpar *= (lambdaval+4.0)*scra; <<
432 } <<
433 // return with the selected GS angular distr <<
434 return dtr; <<
435 } 670 }
436 671
437 672
438 void G4GoudsmitSaundersonTable::LoadMSCData() << 673 std::string pathString(path);
439 gGSMSCAngularDistributions1.resize(gLAMBNUM* << 674 sprintf(basename,"inverseq2CDFII");
440 const G4String str1 = G4EmParameters::Instan << 675 for(int k = 0; k < fgNumLambdas; ++k){
441 for (G4int il=0; il<gLAMBNUM; ++il) { << 676 char fname[512];
442 G4String fname = str1 + std::to_string(il) << 677 sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
443 std::ifstream infile(fname,std::ios::in); << 678 std::ifstream infile(fname,std::ios::in);
444 if (!infile.is_open()) { << 679 if(!infile.is_open()){
445 G4String msgc = "Cannot open file: " + f << 680 char msgc[512];
446 G4Exception("G4GoudsmitSaundersonTable:: << 681 sprintf(msgc,"Cannot open file: %s .",fname);
447 FatalException, msgc.c_str()); << 682 G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
448 return; << 683 FatalException,
449 } << 684 msgc);
450 for (G4int iq=0; iq<gQNUM1; ++iq) { << 685 return;
451 auto gsd = new GSMSCAngularDtr(); << 686 }
452 infile >> gsd->fNumData; <<
453 gsd->fUValues = new G4double[gsd->fNumDa <<
454 gsd->fParamA = new G4double[gsd->fNumDa <<
455 gsd->fParamB = new G4double[gsd->fNumDa <<
456 G4double ddummy; <<
457 infile >> ddummy; infile >> ddummy; <<
458 for (G4int i=0; i<gsd->fNumData; ++i) { <<
459 infile >> gsd->fUValues[i]; <<
460 infile >> gsd->fParamA[i]; <<
461 infile >> gsd->fParamB[i]; <<
462 } <<
463 gGSMSCAngularDistributions1[il*gQNUM1+iq <<
464 } <<
465 infile.close(); <<
466 } <<
467 // <<
468 // second grid <<
469 gGSMSCAngularDistributions2.resize(gLAMBNUM* <<
470 const G4String str2 = G4EmParameters::Instan <<
471 for (G4int il=0; il<gLAMBNUM; ++il) { <<
472 G4String fname = str2 + std::to_string(il) <<
473 std::ifstream infile(fname,std::ios::in); <<
474 if (!infile.is_open()) { <<
475 G4String msgc = "Cannot open file: " + f <<
476 G4Exception("G4GoudsmitSaundersonTable:: <<
477 FatalException, msgc.c_str()); <<
478 return; <<
479 } <<
480 for (G4int iq=0; iq<gQNUM2; ++iq) { <<
481 G4int numData; <<
482 infile >> numData; <<
483 if (numData>1) { <<
484 auto gsd = new GSMSCAngularDtr(); <<
485 gsd->fNumData = numData; <<
486 gsd->fUValues = new G4double[gsd->fNum <<
487 gsd->fParamA = new G4double[gsd->fNum <<
488 gsd->fParamB = new G4double[gsd->fNum <<
489 double ddummy; <<
490 infile >> ddummy; infile >> ddummy; <<
491 for (G4int i=0; i<gsd->fNumData; ++i) <<
492 infile >> gsd->fUValues[i]; <<
493 infile >> gsd->fParamA[i]; <<
494 infile >> gsd->fParamB[i]; <<
495 } <<
496 gGSMSCAngularDistributions2[il*gQNUM2+ <<
497 } else { <<
498 gGSMSCAngularDistributions2[il*gQNUM2+ <<
499 } <<
500 } <<
501 infile.close(); <<
502 } <<
503 } <<
504 687
505 // samples cost in single scattering based on << 688 int limk = k*fgNumUvalues*fgNumLamG1II;
506 // (with Mott-correction if it was requested) << 689 for(int i = 0; i < fgNumUvalues; ++i)
507 G4double G4GoudsmitSaundersonTable::SingleScat << 690 for(int j = 0; j < fgNumLamG1II+1; ++j)
508 << 691 if(j==0) infile >> dum;
509 << 692 else infile >> fgInverseQ2CDFsII[limk+(j-1)*fgNumUvalues+i];
510 G4double rand1 = G4UniformRand(); << 693
511 // sample cost from the Screened-Rutherford << 694 infile.close();
512 G4double cost = 1.-2.0*scra*rand1/(1.0-rand << 695 }
513 // Mott-correction if it was requested by th << 696
514 if (fIsMottCorrection) { << 697
515 static const G4int nlooplim = 1000; // re << 698 sprintf(basename,"inverseq2CDFRatInterpPAII");
516 G4int nloop = 0 ; // loop counter << 699 for(int k = 0; k < fgNumLambdas; ++k){
517 G4int ekindx = -1 ; // evaluate only << 700 char fname[512];
518 G4int deltindx = 0 ; // single scatter << 701 sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
519 G4double q1 = 0.; // not used when << 702 std::ifstream infile(fname,std::ios::in);
520 // computing Mott rejection function value << 703 if(!infile.is_open()){
521 G4double val = fMottCorrection->GetMo << 704 char msgc[512];
522 << 705 sprintf(msgc,"Cannot open file: %s .",fname);
523 while (G4UniformRand()>val && ++nloop<nloo << 706 G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
524 // sampling cos(theta) from the Screened << 707 FatalException,
525 rand1 = G4UniformRand(); << 708 msgc);
526 cost = 1.-2.0*scra*rand1/(1.0-rand1+scr << 709 return;
527 // computing Mott rejection function val << 710 }
528 val = fMottCorrection->GetMottRejectio <<
529 <<
530 }; <<
531 } <<
532 return cost; <<
533 } <<
534 711
>> 712 int limk = k*fgNumUvalues*fgNumLamG1II;
>> 713 for(int i = 0; i < fgNumUvalues; ++i)
>> 714 for(int j = 0; j < fgNumLamG1II+1; ++j)
>> 715 if(j==0) infile >> dum;
>> 716 else infile >> fgInterParamsA2II[limk+(j-1)*fgNumUvalues+i];
>> 717
>> 718 infile.close();
>> 719 }
>> 720
>> 721 sprintf(basename,"inverseq2CDFRatInterpPBII");
>> 722 for(int k = 0; k < fgNumLambdas; ++k){
>> 723 char fname[512];
>> 724 sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
>> 725 std::ifstream infile(fname,std::ios::in);
>> 726 if(!infile.is_open()){
>> 727 char msgc[512];
>> 728 sprintf(msgc,"Cannot open file: %s .",fname);
>> 729 G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
>> 730 FatalException,
>> 731 msgc);
>> 732 return;
>> 733 }
535 734
536 void G4GoudsmitSaundersonTable::GetMottCorrec << 735 int limk = k*fgNumUvalues*fgNumLamG1II;
537 << 736 for(int i = 0; i < fgNumUvalues; ++i)
538 << 737 for(int j = 0; j < fgNumLamG1II+1; ++j)
539 if (fIsMottCorrection) { << 738 if(j==0) infile >> dum;
540 fMottCorrection->GetMottCorrectionFactors( << 739 else infile >> fgInterParamsB2II[limk+(j-1)*fgNumUvalues+i];
541 } <<
542 } <<
543 740
>> 741 infile.close();
>> 742 }
>> 743 }
544 744
>> 745 ////////////////////////////////////////////////////////////////////////////////
545 // compute material dependent Moliere MSC para 746 // compute material dependent Moliere MSC parameters at initialisation
546 void G4GoudsmitSaundersonTable::InitMoliereMSC << 747 void G4GoudsmitSaundersonTable::InitMoliereMSCParams(){
547 const G4double const1 = 7821.6; // [ 748 const G4double const1 = 7821.6; // [cm2/g]
548 const G4double const2 = 0.1569; // [ 749 const G4double const2 = 0.1569; // [cm2 MeV2 / g]
549 const G4double finstrc2 = 5.325135453E-5; / 750 const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
550 751
551 G4MaterialTable* theMaterialTable = G4Mater << 752 G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable();
552 // get number of materials in the table 753 // get number of materials in the table
553 std::size_t numMaterials = theMaterialTable << 754 unsigned int numMaterials = theMaterialTable->size();
554 // make sure that we have long enough vecto 755 // make sure that we have long enough vectors
555 if(gMoliereBc.size()<numMaterials) { << 756 if(!fgMoliereBc) {
556 gMoliereBc.resize(numMaterials); << 757 fgMoliereBc = new std::vector<G4double>(numMaterials);
557 gMoliereXc2.resize(numMaterials); << 758 fgMoliereXc2 = new std::vector<G4double>(numMaterials);
558 } << 759 } else {
559 G4double xi = 1.0; << 760 fgMoliereBc->resize(numMaterials);
560 G4int maxZ = 200; << 761 fgMoliereXc2->resize(numMaterials);
561 if (fIsMottCorrection || fIsPWACorrection) <<
562 // xi = 1.0; <= always set to 1 from n <<
563 maxZ = G4GSMottCorrection::GetMaxZet(); <<
564 } 762 }
565 // << 763
566 for (std::size_t imat=0; imat<numMaterials; << 764 const G4double xims = 0.0; // use xi_MS = 0.0 value now. We will see later on.
567 const G4Material* theMaterial = << 765 for(unsigned int imat = 0; imat < numMaterials; ++imat) {
568 const G4ElementVector* theElemVect = << 766 const G4Material *theMaterial = (*theMaterialTable)[imat];
569 const G4int numelems = << 767 const G4ElementVector *theElemVect = theMaterial->GetElementVector();
570 // << 768
571 const G4double* theNbAtomsPerVolVe << 769 const G4int numelems = theMaterial->GetNumberOfElements();
572 G4double theTotNbAtomsPerVo << 770 const G4double* theFractionVect = theMaterial->GetFractionVector();
573 // <<
574 G4double zs = 0.0; 771 G4double zs = 0.0;
575 G4double zx = 0.0; 772 G4double zx = 0.0;
576 G4double ze = 0.0; 773 G4double ze = 0.0;
577 G4double sa = 0.0; << 774
578 // <<
579 for(G4int ielem = 0; ielem < numelems; ie 775 for(G4int ielem = 0; ielem < numelems; ielem++) {
580 G4double zet = (*theElemVect)[ielem]->G << 776 G4double izet = (*theElemVect)[ielem]->GetZ();
581 if (zet>maxZ) { <<
582 zet = (G4double)maxZ; <<
583 } <<
584 G4double iwa = (*theElemVect)[ielem]-> 777 G4double iwa = (*theElemVect)[ielem]->GetN();
585 G4double ipz = theNbAtomsPerVolVect[ie << 778 // G4int ipz = theAtomsVect[ielem];
586 G4double dum = ipz*zet*(zet+xi); << 779 G4double ipz = theFractionVect[ielem];
587 zs += dum; << 780
588 ze += dum*(-2.0/3.0)*G4Log(ze << 781 G4double dum = ipz/iwa*izet*(izet+xims);
589 zx += dum*G4Log(1.0+3.34*fins << 782 zs += dum;
590 sa += ipz*iwa; << 783 ze += dum*(-2.0/3.0)*G4Log(izet);
>> 784 zx += dum*G4Log(1.0+3.34*finstrc2*izet*izet);
>> 785 // wa += ipz*iwa;
591 } 786 }
592 G4double density = theMaterial->GetDensit 787 G4double density = theMaterial->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
593 // <<
594 gMoliereBc[theMaterial->GetIndex()] = co <<
595 gMoliereXc2[theMaterial->GetIndex()] = co <<
596 // change to Geant4 internal units of 1/l <<
597 gMoliereBc[theMaterial->GetIndex()] *= 1 <<
598 gMoliereXc2[theMaterial->GetIndex()] *= C <<
599 } <<
600 } <<
601 <<
602 788
603 // this method is temporary, will be removed/r << 789 (*fgMoliereBc)[theMaterial->GetIndex()] = const1*density*zs*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
604 G4double G4GoudsmitSaundersonTable::ComputeSca << 790 (*fgMoliereXc2)[theMaterial->GetIndex()] = const2*density*zs; // [MeV2/cm]
605 G4int imc = matcut->GetIndex(); <<
606 G4double corFactor = 1.0; <<
607 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin< <<
608 return corFactor; <<
609 } <<
610 // get the scattering power correction facto <<
611 G4double lekin = G4Log(ekin); <<
612 G4double remaining = (lekin-fSCPCPerMatCuts <<
613 std::size_t lindx = (std::size_t)remaining <<
614 remaining -= lindx; <<
615 std::size_t imax = fSCPCPerMatCuts[imc]-> <<
616 if (lindx>=imax) { <<
617 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i <<
618 } else { <<
619 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l <<
620 } <<
621 return corFactor; <<
622 } <<
623 791
>> 792 // change to Geant4 internal units of 1/length and energ2/length
>> 793 (*fgMoliereBc)[theMaterial->GetIndex()] *= 1.0/CLHEP::cm;
>> 794 (*fgMoliereXc2)[theMaterial->GetIndex()] *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
>> 795 }
624 796
625 void G4GoudsmitSaundersonTable::InitSCPCorrect <<
626 // get the material-cuts table <<
627 G4ProductionCutsTable *thePCTable = G4Produc <<
628 std::size_t numMatCuts = thePCTab <<
629 // clear container if any <<
630 for (std::size_t imc=0; imc<fSCPCPerMatCuts. <<
631 if (fSCPCPerMatCuts[imc]) { <<
632 fSCPCPerMatCuts[imc]->fVSCPC.clear(); <<
633 delete fSCPCPerMatCuts[imc]; <<
634 fSCPCPerMatCuts[imc] = nullptr; <<
635 } <<
636 } <<
637 // <<
638 // set size of the container and create the <<
639 fSCPCPerMatCuts.resize(numMatCuts,nullptr); <<
640 // loop over the material-cuts and create sc <<
641 for (G4int imc=0; imc<(G4int)numMatCuts; ++i <<
642 const G4MaterialCutsCouple *matCut = theP <<
643 // get e- production cut in the current ma <<
644 G4double limit; <<
645 G4double ecut; <<
646 if (fIsElectron) { <<
647 ecut = (*(thePCTable->GetEnergyCutsVect <<
648 limit = 2.*ecut; <<
649 } else { <<
650 ecut = (*(thePCTable->GetEnergyCutsVect <<
651 limit = ecut; <<
652 } <<
653 G4double min = std::max(limit,fLowEnergyLi <<
654 G4double max = fHighEnergyLimit; <<
655 if (min>=max) { <<
656 fSCPCPerMatCuts[imc] = new SCPCorrection <<
657 fSCPCPerMatCuts[imc]->fIsUse = false; <<
658 fSCPCPerMatCuts[imc]->fPrCut = min; <<
659 continue; <<
660 } <<
661 G4int numEbins = fNumSPCEbinPerDec*G <<
662 numEbins = std::max(numEbins,3 <<
663 G4double lmin = G4Log(min); <<
664 G4double ldel = G4Log(max/min)/(num <<
665 fSCPCPerMatCuts[imc] = new SCPCorrection() <<
666 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi <<
667 fSCPCPerMatCuts[imc]->fIsUse = true; <<
668 fSCPCPerMatCuts[imc]->fPrCut = min; <<
669 fSCPCPerMatCuts[imc]->fLEmin = lmin; <<
670 fSCPCPerMatCuts[imc]->fILDel = 1./ldel; <<
671 for (G4int ie=0; ie<numEbins; ++ie) { <<
672 G4double ekin = G4Exp(lmin+ie*ldel); <<
673 G4double scpCorr = 1.0; <<
674 // compute correction factor: I.Kawrakow <<
675 if (ie>0) { <<
676 G4double tau = ekin/CLHEP::electr <<
677 G4double tauCut = ecut/CLHEP::electr <<
678 // Moliere's screening parameter <<
679 G4int matindx = (G4int)matCut->Get <<
680 G4double A = GetMoliereXc2(mati <<
681 G4double gr = (1.+2.*A)*G4Log(1. <<
682 G4double dum0 = (tau+2.)/(tau+1.); <<
683 G4double dum1 = tau+1.; <<
684 G4double gm = G4Log(0.5*tau/tauC <<
685 - 0.25*(tau+2.)*( <<
686 G4Log((tau+4.)*( <<
687 + 0.5*(tau-2*tauCu <<
688 if (gm<gr) { <<
689 gm = gm/gr; <<
690 } else { <<
691 gm = 1.; <<
692 } <<
693 G4double z0 = matCut->GetMaterial()-> <<
694 scpCorr = 1.-gm*z0/(z0*(z0+1.)); <<
695 } <<
696 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo <<
697 } <<
698 } <<
699 } 797 }
700 798