Geant4  10.02
G4GoudsmitSaundersonTable.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // $Id: G4GoudsmitSaundersonTable.cc 93663 2015-10-28 09:50:49Z gcosmo $
27 //
28 // -----------------------------------------------------------------------------
29 //
30 // GEANT4 Class implementation file
31 //
32 // File name: G4GoudsmitSaundersonTable
33 //
34 // Author: Mihaly Novak / (Omrane Kadri)
35 //
36 // Creation date: 20.02.2009
37 //
38 // Class description:
39 // Class to handle multiple scattering angular distributions precomputed by
40 // using Kawrakow-Bielajew Goudsmit-Saunderson MSC model based on the screened
41 // Rutherford DCS for elastic scattering of electrons/positrons [1,2]. This
42 // class is used by G4GoudsmitSaundersonMscModel to sample the angular
43 // deflection of electrons/positrons after travelling a given path.
44 //
45 // Modifications:
46 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
47 // 26.08.2009 O.Kadri: avoiding unuseful calculations and optimizing the root
48 // finding parameter error's within SampleTheta method
49 // 08.02.2010 O.Kadri: reduce delared variables; reduce error of finding root
50 // in secant method
51 // 26.03.2010 O.Kadri: minimum of used arrays in computation within the dichotomie
52 // finding method the error was the lowest value of uvalues
53 // 12.05.2010 O.Kadri: changing of sqrt((b-a)*(b-a)) with fabs(b-a)
54 // 18.05.2015 M. Novak This class has been completely replaced (only the original
55 // class name was kept; class description was also inserted):
56 // A new version of Kawrakow-Bielajew Goudsmit-Saunderson MSC model
57 // based on the screened Rutherford DCS for elastic scattering of
58 // electrons/positrons has been introduced[1,2]. The corresponding MSC
59 // angular distributions over a 2D parameter grid have been recomputed
60 // and the CDFs are now stored in a variable transformed (smooth) form
61 // together with the corresponding rational interpolation parameters.
62 // The new version is several times faster, more robust and accurate
63 // compared to the earlier version (G4GoudsmitSaundersonMscModel class
64 // that use these data has been also completely replaced)
65 //
66 // References:
67 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-208
68 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(1998) 325-336
69 //
70 // -----------------------------------------------------------------------------
71 
72 #include "G4GoudsmitSaundersonTable.hh"
73 
74 #include "Randomize.hh"
75 #include "G4PhysicalConstants.hh"
76 #include "G4MaterialTable.hh"
77 #include "G4Material.hh"
78 #include "G4Log.hh"
79 #include "G4Exp.hh"
80 
81 #include <fstream>
82 #include <cstdlib>
83 #include <cstdio>
84 #include <cmath>
85 
86 
87 const G4double G4GoudsmitSaundersonTable::fgLambdaValues[] ={
88  1.000000000000e+00, 1.165914401180e+00, 1.359356390879e+00, 1.584893192461e+00,
89  1.847849797422e+00, 2.154434690032e+00, 2.511886431510e+00, 2.928644564625e+00,
90  3.414548873834e+00, 3.981071705535e+00, 4.641588833613e+00, 5.411695265465e+00,
91  6.309573444802e+00, 7.356422544596e+00, 8.576958985909e+00, 1.000000000000e+01,
92  1.165914401180e+01, 1.359356390879e+01, 1.584893192461e+01, 1.847849797422e+01,
93  2.154434690032e+01, 2.511886431510e+01, 2.928644564625e+01, 3.414548873834e+01,
94  3.981071705535e+01, 4.641588833613e+01, 5.411695265465e+01, 6.309573444802e+01,
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  };
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 
288 void G4GoudsmitSaundersonTable::Initialise(){
289  if(!fgIsInitialised){
290  LoadMSCData();
291  LoadMSCDataII();
292  fgIsInitialised = TRUE;
293  }
294 
295  InitMoliereMSCParams();
296 }
297 
299 G4GoudsmitSaundersonTable::~G4GoudsmitSaundersonTable(){
300  if(fgMoliereBc){
301  delete fgMoliereBc;
302  fgMoliereBc = 0;
303  }
304  if(fgMoliereXc2){
305  delete fgMoliereXc2;
306  fgMoliereXc2 = 0;
307  }
308  fgIsInitialised = FALSE;
309 }
310 
312 // sample cos(theta) i.e. angular deflection from the precomputed angular
313 // distributions in the real multiple scattering case
314 G4double G4GoudsmitSaundersonTable::SampleCosTheta(G4double lambdavalue, G4double lamG1value,
315  G4double screeningparam, G4double rndm1, G4double rndm2,
316  G4double rndm3){
317  const G4double lnl0 = 0.0; //log(fgLambdaValues[0]);
318  const G4double invlnl = 6.5144172285e+00; //1./log(fgLambdaValues[i+1]/fgLambdaValues[i])
319  G4double logLambdaVal = G4Log(lambdavalue);
320  G4int lambdaindx = (G4int)((logLambdaVal-lnl0)*invlnl);
321 
322  const G4double dd = 1.0/0.0499; // delta
323  G4int lamg1indx = (G4int)((lamG1value-fgLamG1Values[0])*dd);
324 
325  G4double probOfLamJ = G4Log(fgLambdaValues[lambdaindx+1]/lambdavalue)*invlnl;
326  if(rndm1 > probOfLamJ)
327  ++lambdaindx;
328 
329  // store 1.0/(fgLamG1Values[lamg1indx+1]-fgLamG1Values[lamg1indx]) const value
330  const G4double onePerLamG1Diff = dd;
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);
352 
353  // transform back u to cos(theta) :
354  // -compute parameter value a = w2*screening_parameter
355  G4double t = logLambdaVal;
356  G4double dum;
357  if(lambdavalue >= 10.0)
358  dum = 0.5*(-2.77164+t*( 2.94874-t*(0.1535754-t*0.00552888) ));
359  else
360  dum = 0.5*(1.347+t*(0.209364-t*(0.45525-t*(0.50142-t*0.081234))));
361 
362  G4double a = dum*(lambdavalue+4.0)*screeningparam;
363 
364  // this is the sampled cos(theta)
365  return (2.0*a*sample)/(1.0-sample+a);
366 }
367 
369 // sample cos(theta) i.e. angular deflection from the precomputed angular
370 // distributions in the real multiple scattering case (For the second grid)
371 G4double G4GoudsmitSaundersonTable::SampleCosThetaII(G4double lambdavalue, G4double lamG1value,
372  G4double screeningparam, G4double rndm1, G4double rndm2,
373  G4double rndm3){
374  const G4double lnl0 = 0.0; //log(fgLambdaValues[0]);
375  const G4double invlnl = 6.5144172285e+00; //1./log(fgLambdaValues[i+1]/fgLambdaValues[i])
376  G4double logLambdaVal = G4Log(lambdavalue);
377  G4int lambdaindx = (G4int)((logLambdaVal-lnl0)*invlnl);
378 
379  const G4double dd = 1.0/0.333; // delta
380  G4int lamg1indx = (G4int)((lamG1value-fgLamG1ValuesII[0])*dd);
381 
382  G4double probOfLamJ = G4Log(fgLambdaValues[lambdaindx+1]/lambdavalue)*invlnl;
383  if(rndm1 > probOfLamJ)
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  G4int begin = lambdaindx*(fgNumLamG1II*fgNumUvalues)+lamg1indx*fgNumUvalues;
393 
394  // sample bin
395  G4int indx = (G4int)(rndm3*100); // value/delatU ahol deltaU = 0.01; find u-indx
396  G4double cp1 = fgUValues[indx];
397 // G4double cp2 = fgUValues[indx+1];
398  indx +=begin;
399 
400  G4double u1 = fgInverseQ2CDFsII[indx];
401  G4double u2 = fgInverseQ2CDFsII[indx+1];
402 
403  G4double dum1 = rndm3 - cp1;
404  const G4double dum2 = 0.01;//cp2-cp1; // this is constans 0.01
405  const G4double dum22 = 0.0001;
406  // rational interpolation of the inverse CDF
407  G4double sample= u1+(1.0+fgInterParamsA2II[indx]+fgInterParamsB2II[indx])*dum1*dum2/(dum22+
408  dum1*dum2*fgInterParamsA2II[indx]+fgInterParamsB2II[indx]*dum1*dum1)*(u2-u1);
409 
410  // transform back u to cos(theta) :
411  // -compute parameter value a = w2*screening_parameter
412  G4double t = logLambdaVal;
413  G4double dum = 0.;
414  if(lambdavalue >= 10.0)
415  dum = 0.5*(-2.77164+t*( 2.94874-t*(0.1535754-t*0.00552888) ));
416  else
417  dum = 0.5*(1.347+t*(0.209364-t*(0.45525-t*(0.50142-t*0.081234))));
418 
419  G4double a = dum*(lambdavalue+4.0)*screeningparam;
420 
421  // this is the sampled cos(theta)
422  return (2.0*a*sample)/(1.0-sample+a);
423 }
424 
426 // samplig multiple scattering angles cos(theta) and sin(thata)
427 // including no-scattering, single, "few" scattering cases as well
428 void G4GoudsmitSaundersonTable::Sampling(G4double lambdavalue, G4double lamG1value,
429  G4double scrPar, G4double &cost, G4double &sint){
430  G4double rand0 = G4UniformRand();
431  G4double expn = G4Exp(-lambdavalue);
432 
433  //
434  // no scattering case
435  //
436  if(rand0 < expn) {
437  cost = 1.0;
438  sint = 0.0;
439  return;
440  }
441 
442  //
443  // single scattering case : sample (analitically) from the single scattering PDF
444  // that corresponds to the modified-Wentzel DCS i.e.
445  // Wentzel DCS that gives back the PWA first-transport mfp
446  //
447  if( rand0 < (1.+lambdavalue)*expn) {
448  G4double rand1 = G4UniformRand();
449  // sampling 1-cos(theta)
450  G4double dum0 = 2.0*scrPar*rand1/(1.0 - rand1 + scrPar);
451  // add protections
452  if(dum0 < 0.0) dum0 = 0.0;
453  else if(dum0 > 2.0 ) dum0 = 2.0;
454  // compute cos(theta) and sin(theta)
455  cost = 1.0 - dum0;
456  sint = std::sqrt(dum0*(2.0 - dum0));
457  return;
458  }
459 
460  //
461  // handle this case:
462  // -lambdavalue < 1 i.e. mean #elastic events along the step is < 1 but
463  // the currently sampled case is not 0 or 1 scattering. [Our minimal
464  // lambdavalue (that we have precomputed, transformed angular distributions
465  // stored in a form of equally probabe intervalls together with rational
466  // interp. parameters) is 1.]
467  // -probability of having n elastic events follows Poisson stat. with
468  // lambdavalue parameter.
469  // -the max. probability (when lambdavalue=1) of having more than one
470  // elastic events is 0.2642411 and the prob of having 2,3,..,n elastic
471  // events decays rapidly with n. So set a max n to 10.
472  // -sampling of this cases is done in a one-by-one single elastic event way
473  // where the current #elastic event is sampled from the Poisson distr.
474  // -(other possibility would be to use lambdavalue=1 distribution because
475  // they are very close)
476  if(lambdavalue < 1.0) {
477  G4double prob, cumprob;
478  prob = cumprob = expn;
479  G4double curcost,cursint;
480  // init cos(theta) and sin(theta) to the zero scattering values
481  cost = 1.0;
482  sint = 0.0;
483 
484  for(G4int iel = 1; iel < 10; ++iel){
485  // prob of having iel scattering from Poisson
486  prob *= lambdavalue/(G4double)iel;
487  cumprob += prob;
488  //
489  //sample cos(theta) from the singe scattering pdf
490  //
491  G4double rand1 = G4UniformRand();
492  // sampling 1-cos(theta)
493  G4double dum0 = 2.0*scrPar*rand1/(1.0 - rand1 + scrPar);
494  // compute cos(theta) and sin(theta)
495  curcost = 1.0 - dum0;
496  cursint = dum0*(2.0 - dum0);
497  //
498  // if we got current deflection that is not too small
499  // then update cos(theta) sin(theta)
500  if(cursint > 1.0e-20){
501  cursint = std::sqrt(cursint);
502  G4double curphi = CLHEP::twopi*G4UniformRand();
503  cost = cost*curcost - sint*cursint*std::cos(curphi);
504  sint = std::sqrt(std::max(0.0, (1.0-cost)*(1.0+cost)));
505  }
506  //
507  // check if we have done enough scattering i.e. sampling from the Poisson
508  if( rand0 < cumprob)
509  return;
510  }
511  // if reached the max iter i.e. 10
512  return;
513  }
514 
515 
517  // IT IS DONE NOW IN THE MODEL
518  //if(lamG1value>fgLamG1Values[fgNumLamG1-1]) {
519  // cost = 1.-2.*G4UniformRand();
520  // sint = std::sqrt((1.-cost)*(1.+cost));
521  // return;
522  //}
523 
524 
525  //
526  // multiple scattering case with lambdavalue >= 1:
527  // - use the precomputed and transformed Goudsmit-Saunderson angular
528  // distributions that are stored in a form of equally probabe intervalls
529  // together with the corresponding rational interp. parameters
530  //
531  // add protections
532  if(lamG1value<fgLamG1Values[0]) lamG1value = fgLamG1Values[0];
533  if(lamG1value>fgLamG1ValuesII[fgNumLamG1II-1]) lamG1value = fgLamG1ValuesII[fgNumLamG1II-1];
534 
535  // sample 1-cos(theta) :: depending on the grids
536  G4double dum0 = 0.0;
537  if(lamG1value<fgLamG1Values[fgNumLamG1-1]) {
538  dum0 = SampleCosTheta(lambdavalue, lamG1value, scrPar, G4UniformRand(),
539  G4UniformRand(), G4UniformRand());
540 
541  } else {
542  dum0 = SampleCosThetaII(lambdavalue, lamG1value, scrPar, G4UniformRand(),
543  G4UniformRand(), G4UniformRand());
544  }
545 
546  // add protections
547  if(dum0 < 0.0) dum0 = 0.0;
548  if(dum0 > 2.0 ) dum0 = 2.0;
549 
550  // compute cos(theta) and sin(theta)
551  cost = 1.0-dum0;
552  sint = std::sqrt(dum0*(2.0 - dum0));
553 }
554 
556 G4double G4GoudsmitSaundersonTable::GetScreeningParam(G4double G1){
557  if(G1<fgG1Values[0]) return fgScreeningParam[0];
558  if(G1>fgG1Values[fgNumScreeningParams-1]) return fgScreeningParam[fgNumScreeningParams-1];
559  // normal case
560  const G4double lng10 = -2.30258509299405e+01; //log(fgG1Values[0]);
561  const G4double invlng1 = 6.948711710452e+00; //1./log(fgG1Values[i+1]/fgG1Values[i])
562  G4double logG1 = G4Log(G1);
563  G4int k = (G4int)((logG1-lng10)*invlng1);
564  return G4Exp(logG1*fgSrcAValues[k])*fgSrcBValues[k]; // log-log linear
565 }
566 
568 void G4GoudsmitSaundersonTable::LoadMSCData(){
569  G4double dum;
570  char basename[512];
571  char* path = getenv("G4LEDATA");
572  if (!path) {
573  G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
574  FatalException,
575  "Environment variable G4LEDATA not defined");
576  return;
577  }
578 
579  std::string pathString(path);
580  sprintf(basename,"inverseq2CDF");
581  for(int k = 0; k < fgNumLambdas; ++k){
582  char fname[512];
583  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
584  std::ifstream infile(fname,std::ios::in);
585  if(!infile.is_open()){
586  char msgc[512];
587  sprintf(msgc,"Cannot open file: %s .",fname);
588  G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
589  FatalException,
590  msgc);
591  return;
592  }
593 
594  int limk = k*fgNumUvalues*fgNumLamG1;
595  for(int i = 0; i < fgNumUvalues; ++i)
596  for(int j = 0; j < fgNumLamG1+1; ++j)
597  if(j==0) infile >> dum;
598  else infile >> fgInverseQ2CDFs[limk+(j-1)*fgNumUvalues+i];
599 
600  infile.close();
601  }
602 
603 
604  sprintf(basename,"inverseq2CDFRatInterpPA");
605  for(int k = 0; k < fgNumLambdas; ++k){
606  char fname[512];
607  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
608  std::ifstream infile(fname,std::ios::in);
609  if(!infile.is_open()){
610  char msgc[512];
611  sprintf(msgc,"Cannot open file: %s .",fname);
612  G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
613  FatalException,
614  msgc);
615  return;
616  }
617 
618  int limk = k*fgNumUvalues*fgNumLamG1;
619  for(int i = 0; i < fgNumUvalues; ++i)
620  for(int j = 0; j < fgNumLamG1+1; ++j)
621  if(j==0) infile >> dum;
622  else infile >> fgInterParamsA2[limk+(j-1)*fgNumUvalues+i];
623 
624  infile.close();
625  }
626 
627  sprintf(basename,"inverseq2CDFRatInterpPB");
628  for(int k = 0; k < fgNumLambdas; ++k){
629  char fname[512];
630  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
631  std::ifstream infile(fname,std::ios::in);
632  if(!infile.is_open()){
633  char msgc[512];
634  sprintf(msgc,"Cannot open file: %s .",fname);
635  G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
636  FatalException,
637  msgc);
638  return;
639  }
640 
641  int limk = k*fgNumUvalues*fgNumLamG1;
642  for(int i = 0; i < fgNumUvalues; ++i)
643  for(int j = 0; j < fgNumLamG1+1; ++j)
644  if(j==0) infile >> dum;
645  else infile >> fgInterParamsB2[limk+(j-1)*fgNumUvalues+i];
646 
647  infile.close();
648  }
649 }
650 
652 // Load the second grid data
654 void G4GoudsmitSaundersonTable::LoadMSCDataII(){
655 G4double dum;
656 char basename[512];
657 char* path = getenv("G4LEDATA");
658 if (!path) {
659  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
660  FatalException,
661  "Environment variable G4LEDATA not defined");
662  return;
663 }
664 
665 
666  std::string pathString(path);
667  sprintf(basename,"inverseq2CDFII");
668  for(int k = 0; k < fgNumLambdas; ++k){
669  char fname[512];
670  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
671  std::ifstream infile(fname,std::ios::in);
672  if(!infile.is_open()){
673  char msgc[512];
674  sprintf(msgc,"Cannot open file: %s .",fname);
675  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
676  FatalException,
677  msgc);
678  return;
679  }
680 
681  int limk = k*fgNumUvalues*fgNumLamG1II;
682  for(int i = 0; i < fgNumUvalues; ++i)
683  for(int j = 0; j < fgNumLamG1II+1; ++j)
684  if(j==0) infile >> dum;
685  else infile >> fgInverseQ2CDFsII[limk+(j-1)*fgNumUvalues+i];
686 
687  infile.close();
688  }
689 
690 
691  sprintf(basename,"inverseq2CDFRatInterpPAII");
692  for(int k = 0; k < fgNumLambdas; ++k){
693  char fname[512];
694  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
695  std::ifstream infile(fname,std::ios::in);
696  if(!infile.is_open()){
697  char msgc[512];
698  sprintf(msgc,"Cannot open file: %s .",fname);
699  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
700  FatalException,
701  msgc);
702  return;
703  }
704 
705  int limk = k*fgNumUvalues*fgNumLamG1II;
706  for(int i = 0; i < fgNumUvalues; ++i)
707  for(int j = 0; j < fgNumLamG1II+1; ++j)
708  if(j==0) infile >> dum;
709  else infile >> fgInterParamsA2II[limk+(j-1)*fgNumUvalues+i];
710 
711  infile.close();
712  }
713 
714  sprintf(basename,"inverseq2CDFRatInterpPBII");
715  for(int k = 0; k < fgNumLambdas; ++k){
716  char fname[512];
717  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
718  std::ifstream infile(fname,std::ios::in);
719  if(!infile.is_open()){
720  char msgc[512];
721  sprintf(msgc,"Cannot open file: %s .",fname);
722  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
723  FatalException,
724  msgc);
725  return;
726  }
727 
728  int limk = k*fgNumUvalues*fgNumLamG1II;
729  for(int i = 0; i < fgNumUvalues; ++i)
730  for(int j = 0; j < fgNumLamG1II+1; ++j)
731  if(j==0) infile >> dum;
732  else infile >> fgInterParamsB2II[limk+(j-1)*fgNumUvalues+i];
733 
734  infile.close();
735  }
736 }
737 
739 // compute material dependent Moliere MSC parameters at initialisation
740 void G4GoudsmitSaundersonTable::InitMoliereMSCParams(){
741  const G4double const1 = 7821.6; // [cm2/g]
742  const G4double const2 = 0.1569; // [cm2 MeV2 / g]
743  const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
744 
745  G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable();
746  // get number of materials in the table
747  unsigned int numMaterials = theMaterialTable->size();
748  // make sure that we have long enough vectors
749  if(!fgMoliereBc) {
750  fgMoliereBc = new std::vector<G4double>(numMaterials);
751  fgMoliereXc2 = new std::vector<G4double>(numMaterials);
752  } else {
753  fgMoliereBc->resize(numMaterials);
754  fgMoliereXc2->resize(numMaterials);
755  }
756 
757  const G4double xims = 0.0; // use xi_MS = 0.0 value now. We will see later on.
758  for(unsigned int imat = 0; imat < numMaterials; ++imat) {
759  const G4Material *theMaterial = (*theMaterialTable)[imat];
760  const G4ElementVector *theElemVect = theMaterial->GetElementVector();
761 
762  const G4int numelems = theMaterial->GetNumberOfElements();
763  const G4double* theFractionVect = theMaterial->GetFractionVector();
764  G4double zs = 0.0;
765  G4double zx = 0.0;
766  G4double ze = 0.0;
767 
768  for(G4int ielem = 0; ielem < numelems; ielem++) {
769  G4double izet = (*theElemVect)[ielem]->GetZ();
770  G4double iwa = (*theElemVect)[ielem]->GetN();
771 // G4int ipz = theAtomsVect[ielem];
772  G4double ipz = theFractionVect[ielem];
773 
774  G4double dum = ipz/iwa*izet*(izet+xims);
775  zs += dum;
776  ze += dum*(-2.0/3.0)*G4Log(izet);
777  zx += dum*G4Log(1.0+3.34*finstrc2*izet*izet);
778 // wa += ipz*iwa;
779  }
780  G4double density = theMaterial->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
781 
782  (*fgMoliereBc)[theMaterial->GetIndex()] = const1*density*zs*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
783  (*fgMoliereXc2)[theMaterial->GetIndex()] = const2*density*zs; // [MeV2/cm]
784 
785  // change to Geant4 internal units of 1/length and energ2/length
786  (*fgMoliereBc)[theMaterial->GetIndex()] *= 1.0/CLHEP::cm;
787  (*fgMoliereXc2)[theMaterial->GetIndex()] *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
788  }
789 
790 }
G4Exp.hh
cm
static const double cm
Definition: G4SIunits.hh:118
G4GoudsmitSaundersonTable::fgNumLamG1II
static const G4int fgNumLamG1II
number of s/{e}G_{1} $-values
Definition: G4GoudsmitSaundersonTable.hh:127
G4GoudsmitSaundersonTable::fgLambdaValues
static const G4double fgLambdaValues[]
number of s/{e} $-values; size = fgNumLambdas = 76
Definition: G4GoudsmitSaundersonTable.hh:135
MeV
static const double MeV
Definition: G4SIunits.hh:211
G4ElementVector
std::vector< G4Element * > G4ElementVector
Definition: G4ElementVector.hh:44
G4Material::GetIndex
size_t GetIndex() const
Definition: G4Material.hh:262
G4Material::GetMaterialTable
static G4MaterialTable * GetMaterialTable()
Definition: G4Material.cc:589
G4MaterialTable
std::vector< G4Material * > G4MaterialTable
Definition: G4MaterialTable.hh:41
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:180
a
G4double a
Definition: TRTMaterials.hh:39
G4Material::GetElementVector
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
G4int
int G4int
Definition: G4Types.hh:78
G4GoudsmitSaundersonTable::fgNumLambdas
static const G4int fgNumLambdas
size of grids of some parameters
Definition: G4GoudsmitSaundersonTable.hh:125
density
G4double density
Definition: TRTMaterials.hh:39
G4GoudsmitSaundersonTable::SampleCosThetaII
G4double SampleCosThetaII(G4double, G4double, G4double, G4double, G4double, G4double)
Definition: G4GoudsmitSaundersonTable.cc:371
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4bool
bool G4bool
Definition: G4Types.hh:79
G4Material.hh
FALSE
#define FALSE
Definition: globals.hh:52
twopi
static const double twopi
Definition: G4SIunits.hh:75
cm3
static const double cm3
Definition: G4SIunits.hh:120
G4GoudsmitSaundersonTable.hh
TRUE
#define TRUE
Definition: globals.hh:55
Randomize.hh
G4GoudsmitSaundersonTable::Sampling
void Sampling(G4double, G4double, G4double, G4double &, G4double &)
Definition: G4GoudsmitSaundersonTable.cc:428
G4Log.hh
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4PhysicalConstants.hh
G4GoudsmitSaundersonTable::fgMoliereBc
static std::vector< G4double > * fgMoliereBc
the grid of b_lambda_{c} $ and {e} $) under the screened Rutherford cross section approximation...
Definition: G4GoudsmitSaundersonTable.hh:173
G4Log
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
g
static const double g
Definition: G4SIunits.hh:180
FatalException
Definition: G4ExceptionSeverity.hh:60
G4GoudsmitSaundersonTable::fgMoliereXc2
static std::vector< G4double > * fgMoliereXc2
Definition: G4GoudsmitSaundersonTable.hh:174
G4MaterialTable.hh
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
G4double
double G4double
Definition: G4Types.hh:76
G4GoudsmitSaundersonTable::SampleCosTheta
G4double SampleCosTheta(G4double, G4double, G4double, G4double, G4double, G4double)
Definition: G4GoudsmitSaundersonTable.cc:314
G4Material::GetFractionVector
const G4double * GetFractionVector() const
Definition: G4Material.hh:194