Geant4  10.02.p02
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 95473 2016-02-12 09:36:21Z 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  // 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))));
425 
426  G4double a = dum*(lambdavalue+4.0)*screeningparam;
427 
428  // this is the sampled cos(theta)
429  return (2.0*a*sample)/(1.0-sample+a);
430 }
431 
433 // samplig multiple scattering angles cos(theta) and sin(thata)
434 // including no-scattering, single, "few" scattering cases as well
435 void G4GoudsmitSaundersonTable::Sampling(G4double lambdavalue, G4double lamG1value,
436  G4double scrPar, G4double &cost, G4double &sint){
437  G4double rand0 = G4UniformRand();
438  G4double expn = G4Exp(-lambdavalue);
439 
440  //
441  // no scattering case
442  //
443  if(rand0 < expn) {
444  cost = 1.0;
445  sint = 0.0;
446  return;
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
453  //
454  if( rand0 < (1.+lambdavalue)*expn) {
455  G4double rand1 = G4UniformRand();
456  // sampling 1-cos(theta)
457  G4double dum0 = 2.0*scrPar*rand1/(1.0 - rand1 + scrPar);
458  // add protections
459  if(dum0 < 0.0) dum0 = 0.0;
460  else if(dum0 > 2.0 ) dum0 = 2.0;
461  // compute cos(theta) and sin(theta)
462  cost = 1.0 - dum0;
463  sint = std::sqrt(dum0*(2.0 - dum0));
464  return;
465  }
466 
467  //
468  // handle this case:
469  // -lambdavalue < 1 i.e. mean #elastic events along the step is < 1 but
470  // the currently sampled case is not 0 or 1 scattering. [Our minimal
471  // lambdavalue (that we have precomputed, transformed angular distributions
472  // stored in a form of equally probabe intervalls together with rational
473  // interp. parameters) is 1.]
474  // -probability of having n elastic events follows Poisson stat. with
475  // lambdavalue parameter.
476  // -the max. probability (when lambdavalue=1) of having more than one
477  // elastic events is 0.2642411 and the prob of having 2,3,..,n elastic
478  // events decays rapidly with n. So set a max n to 10.
479  // -sampling of this cases is done in a one-by-one single elastic event way
480  // where the current #elastic event is sampled from the Poisson distr.
481  // -(other possibility would be to use lambdavalue=1 distribution because
482  // they are very close)
483  if(lambdavalue < 1.0) {
484  G4double prob, cumprob;
485  prob = cumprob = expn;
486  G4double curcost,cursint;
487  // init cos(theta) and sin(theta) to the zero scattering values
488  cost = 1.0;
489  sint = 0.0;
490 
491  for(G4int iel = 1; iel < 10; ++iel){
492  // prob of having iel scattering from Poisson
493  prob *= lambdavalue/(G4double)iel;
494  cumprob += prob;
495  //
496  //sample cos(theta) from the singe scattering pdf
497  //
498  G4double rand1 = G4UniformRand();
499  // sampling 1-cos(theta)
500  G4double dum0 = 2.0*scrPar*rand1/(1.0 - rand1 + scrPar);
501  // compute cos(theta) and sin(theta)
502  curcost = 1.0 - dum0;
503  cursint = dum0*(2.0 - dum0);
504  //
505  // if we got current deflection that is not too small
506  // then update cos(theta) sin(theta)
507  if(cursint > 1.0e-20){
508  cursint = std::sqrt(cursint);
509  G4double curphi = CLHEP::twopi*G4UniformRand();
510  cost = cost*curcost - sint*cursint*std::cos(curphi);
511  sint = std::sqrt(std::max(0.0, (1.0-cost)*(1.0+cost)));
512  }
513  //
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;
520  }
521 
522 
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 
532  //
533  // multiple scattering case with lambdavalue >= 1:
534  // - use the precomputed and transformed Goudsmit-Saunderson angular
535  // distributions that are stored in a form of equally probabe intervalls
536  // together with the corresponding rational interp. parameters
537  //
538  // add protections
539  if(lamG1value<fgLamG1Values[0]) lamG1value = fgLamG1Values[0];
540  if(lamG1value>fgLamG1ValuesII[fgNumLamG1II-1]) lamG1value = fgLamG1ValuesII[fgNumLamG1II-1];
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());
547 
548  } else {
549  dum0 = SampleCosThetaII(lambdavalue, lamG1value, scrPar, G4UniformRand(),
550  G4UniformRand(), G4UniformRand());
551  }
552 
553  // add protections
554  if(dum0 < 0.0) dum0 = 0.0;
555  if(dum0 > 2.0 ) dum0 = 2.0;
556 
557  // compute cos(theta) and sin(theta)
558  cost = 1.0-dum0;
559  sint = std::sqrt(dum0*(2.0 - dum0));
560 }
561 
563 G4double G4GoudsmitSaundersonTable::GetScreeningParam(G4double G1){
564  if(G1<fgG1Values[0]) return fgScreeningParam[0];
565  if(G1>fgG1Values[fgNumScreeningParams-1]) return fgScreeningParam[fgNumScreeningParams-1];
566  // normal case
567  const G4double lng10 = -2.30258509299405e+01; //log(fgG1Values[0]);
568  const G4double invlng1 = 6.948711710452e+00; //1./log(fgG1Values[i+1]/fgG1Values[i])
569  G4double logG1 = G4Log(G1);
570  G4int k = (G4int)((logG1-lng10)*invlng1);
571  return G4Exp(logG1*fgSrcAValues[k])*fgSrcBValues[k]; // log-log linear
572 }
573 
575 void G4GoudsmitSaundersonTable::LoadMSCData(){
576  G4double dum;
577  char basename[512];
578  char* path = getenv("G4LEDATA");
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  }
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  }
624 
625  int limk = k*fgNumUvalues*fgNumLamG1;
626  for(int i = 0; i < fgNumUvalues; ++i)
627  for(int j = 0; j < fgNumLamG1+1; ++j)
628  if(j==0) infile >> dum;
629  else infile >> fgInterParamsA2[limk+(j-1)*fgNumUvalues+i];
630 
631  infile.close();
632  }
633 
634  sprintf(basename,"inverseq2CDFRatInterpPB");
635  for(int k = 0; k < fgNumLambdas; ++k){
636  char fname[512];
637  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
638  std::ifstream infile(fname,std::ios::in);
639  if(!infile.is_open()){
640  char msgc[512];
641  sprintf(msgc,"Cannot open file: %s .",fname);
642  G4Exception("G4GoudsmitSaundersonTable::LoadMSCData()","em0006",
643  FatalException,
644  msgc);
645  return;
646  }
647 
648  int limk = k*fgNumUvalues*fgNumLamG1;
649  for(int i = 0; i < fgNumUvalues; ++i)
650  for(int j = 0; j < fgNumLamG1+1; ++j)
651  if(j==0) infile >> dum;
652  else infile >> fgInterParamsB2[limk+(j-1)*fgNumUvalues+i];
653 
654  infile.close();
655  }
656 }
657 
659 // Load the second grid data
661 void G4GoudsmitSaundersonTable::LoadMSCDataII(){
662 G4double dum;
663 char basename[512];
664 char* path = getenv("G4LEDATA");
665 if (!path) {
666  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
667  FatalException,
668  "Environment variable G4LEDATA not defined");
669  return;
670 }
671 
672 
673  std::string pathString(path);
674  sprintf(basename,"inverseq2CDFII");
675  for(int k = 0; k < fgNumLambdas; ++k){
676  char fname[512];
677  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
678  std::ifstream infile(fname,std::ios::in);
679  if(!infile.is_open()){
680  char msgc[512];
681  sprintf(msgc,"Cannot open file: %s .",fname);
682  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
683  FatalException,
684  msgc);
685  return;
686  }
687 
688  int limk = k*fgNumUvalues*fgNumLamG1II;
689  for(int i = 0; i < fgNumUvalues; ++i)
690  for(int j = 0; j < fgNumLamG1II+1; ++j)
691  if(j==0) infile >> dum;
692  else infile >> fgInverseQ2CDFsII[limk+(j-1)*fgNumUvalues+i];
693 
694  infile.close();
695  }
696 
697 
698  sprintf(basename,"inverseq2CDFRatInterpPAII");
699  for(int k = 0; k < fgNumLambdas; ++k){
700  char fname[512];
701  sprintf(fname,"%s/msc_GS/%s_%d",path,basename,k);
702  std::ifstream infile(fname,std::ios::in);
703  if(!infile.is_open()){
704  char msgc[512];
705  sprintf(msgc,"Cannot open file: %s .",fname);
706  G4Exception("G4GoudsmitSaundersonTable::LoadMSCDataII()","em0006",
707  FatalException,
708  msgc);
709  return;
710  }
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  }
734 
735  int limk = k*fgNumUvalues*fgNumLamG1II;
736  for(int i = 0; i < fgNumUvalues; ++i)
737  for(int j = 0; j < fgNumLamG1II+1; ++j)
738  if(j==0) infile >> dum;
739  else infile >> fgInterParamsB2II[limk+(j-1)*fgNumUvalues+i];
740 
741  infile.close();
742  }
743 }
744 
746 // compute material dependent Moliere MSC parameters at initialisation
747 void G4GoudsmitSaundersonTable::InitMoliereMSCParams(){
748  const G4double const1 = 7821.6; // [cm2/g]
749  const G4double const2 = 0.1569; // [cm2 MeV2 / g]
750  const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
751 
752  G4MaterialTable *theMaterialTable = G4Material::GetMaterialTable();
753  // get number of materials in the table
754  unsigned int numMaterials = theMaterialTable->size();
755  // make sure that we have long enough vectors
756  if(!fgMoliereBc) {
757  fgMoliereBc = new std::vector<G4double>(numMaterials);
758  fgMoliereXc2 = new std::vector<G4double>(numMaterials);
759  } else {
760  fgMoliereBc->resize(numMaterials);
761  fgMoliereXc2->resize(numMaterials);
762  }
763 
764  const G4double xims = 0.0; // use xi_MS = 0.0 value now. We will see later on.
765  for(unsigned int imat = 0; imat < numMaterials; ++imat) {
766  const G4Material *theMaterial = (*theMaterialTable)[imat];
767  const G4ElementVector *theElemVect = theMaterial->GetElementVector();
768 
769  const G4int numelems = theMaterial->GetNumberOfElements();
770  const G4double* theFractionVect = theMaterial->GetFractionVector();
771  G4double zs = 0.0;
772  G4double zx = 0.0;
773  G4double ze = 0.0;
774 
775  for(G4int ielem = 0; ielem < numelems; ielem++) {
776  G4double izet = (*theElemVect)[ielem]->GetZ();
777  G4double iwa = (*theElemVect)[ielem]->GetN();
778 // G4int ipz = theAtomsVect[ielem];
779  G4double ipz = theFractionVect[ielem];
780 
781  G4double dum = ipz/iwa*izet*(izet+xims);
782  zs += dum;
783  ze += dum*(-2.0/3.0)*G4Log(izet);
784  zx += dum*G4Log(1.0+3.34*finstrc2*izet*izet);
785 // wa += ipz*iwa;
786  }
787  G4double density = theMaterial->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
788 
789  (*fgMoliereBc)[theMaterial->GetIndex()] = const1*density*zs*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
790  (*fgMoliereXc2)[theMaterial->GetIndex()] = const2*density*zs; // [MeV2/cm]
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  }
796 
797 }
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:435
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