| 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 // 26 //
27 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class header file 29 // GEANT4 Class header file
30 // 30 //
31 // 31 //
32 // File name: G4eDPWAElasticDCS 32 // File name: G4eDPWAElasticDCS
33 // 33 //
34 // Author: Mihaly Novak 34 // Author: Mihaly Novak
35 // 35 //
36 // Creation date: 02.07.2020 36 // Creation date: 02.07.2020
37 // 37 //
38 // Modifications: 38 // Modifications:
39 // 39 //
40 // 40 //
41 // ------------------------------------------- 41 // -------------------------------------------------------------------
42 42
43 #include "G4eDPWAElasticDCS.hh" 43 #include "G4eDPWAElasticDCS.hh"
44 #include "G4EmParameters.hh" 44 #include "G4EmParameters.hh"
45 #include "G4Physics2DVector.hh" 45 #include "G4Physics2DVector.hh"
46 46
47 #include "zlib.h" 47 #include "zlib.h"
48 48
49 // 49 //
50 // Global variables: 50 // Global variables:
51 // 51 //
52 G4bool G4eDPWAElasticDCS::gIsGr 52 G4bool G4eDPWAElasticDCS::gIsGridLoaded = false;
53 G4String G4eDPWAElasticDCS::gData 53 G4String G4eDPWAElasticDCS::gDataDirectory = "";
54 // final values of these variables will be set 54 // final values of these variables will be set in LoadGrid() called by master
55 std::size_t G4eDPWAElasticDCS::gNumE 55 std::size_t G4eDPWAElasticDCS::gNumEnergies = 106;
56 std::size_t G4eDPWAElasticDCS::gIndx 56 std::size_t G4eDPWAElasticDCS::gIndxEnergyLim = 35;
57 std::size_t G4eDPWAElasticDCS::gNumT 57 std::size_t G4eDPWAElasticDCS::gNumThetas1 = 247;
58 std::size_t G4eDPWAElasticDCS::gNumT 58 std::size_t G4eDPWAElasticDCS::gNumThetas2 = 128;
59 G4double G4eDPWAElasticDCS::gLogM 59 G4double G4eDPWAElasticDCS::gLogMinEkin = 1.0;
60 G4double G4eDPWAElasticDCS::gInvD 60 G4double G4eDPWAElasticDCS::gInvDelLogEkin = 1.0;
61 // containers for grids: Ekin, mu(t)=0.5[1-cos 61 // containers for grids: Ekin, mu(t)=0.5[1-cos(t)] and u(mu,A)=(A+1)mu/(mu+A)
62 std::vector<G4double> G4eDPWAElasticDCS::gTheE 62 std::vector<G4double> G4eDPWAElasticDCS::gTheEnergies(G4eDPWAElasticDCS::gNumEnergies);
63 std::vector<G4double> G4eDPWAElasticDCS::gTheM 63 std::vector<G4double> G4eDPWAElasticDCS::gTheMus1(G4eDPWAElasticDCS::gNumThetas1);
64 std::vector<G4double> G4eDPWAElasticDCS::gTheM 64 std::vector<G4double> G4eDPWAElasticDCS::gTheMus2(G4eDPWAElasticDCS::gNumThetas2);
65 std::vector<G4double> G4eDPWAElasticDCS::gTheU 65 std::vector<G4double> G4eDPWAElasticDCS::gTheU1(G4eDPWAElasticDCS::gNumThetas1);
66 std::vector<G4double> G4eDPWAElasticDCS::gTheU 66 std::vector<G4double> G4eDPWAElasticDCS::gTheU2(G4eDPWAElasticDCS::gNumThetas2);
67 // abscissas and weights of an 8 point Gauss-L 67 // abscissas and weights of an 8 point Gauss-Legendre quadrature
68 // for numerical integration on [0,1] 68 // for numerical integration on [0,1]
69 const G4double G4eDPWAElasticDCS::gXGL[ 69 const G4double G4eDPWAElasticDCS::gXGL[] = {
70 1.98550718E-02, 1.01666761E-01, 2.37233795E- 70 1.98550718E-02, 1.01666761E-01, 2.37233795E-01, 4.08282679E-01,
71 5.91717321E-01, 7.62766205E-01, 8.98333239E- 71 5.91717321E-01, 7.62766205E-01, 8.98333239E-01, 9.80144928E-01
72 }; 72 };
73 const G4double G4eDPWAElasticDCS::gWGL[ 73 const G4double G4eDPWAElasticDCS::gWGL[] = {
74 5.06142681E-02, 1.11190517E-01, 1.56853323E- 74 5.06142681E-02, 1.11190517E-01, 1.56853323E-01, 1.81341892E-01,
75 1.81341892E-01, 1.56853323E-01, 1.11190517E- 75 1.81341892E-01, 1.56853323E-01, 1.11190517E-01, 5.06142681E-02
76 }; 76 };
77 77
78 78
79 // - iselectron : data for e- (for e+ otherw 79 // - iselectron : data for e- (for e+ otherwise)
80 // - isrestricted : sampling of angular deflec 80 // - isrestricted : sampling of angular deflection on restricted interavl is
81 // required (i.e. in case of 81 // required (i.e. in case of mixed-simulation models)
82 G4eDPWAElasticDCS::G4eDPWAElasticDCS(G4bool is 82 G4eDPWAElasticDCS::G4eDPWAElasticDCS(G4bool iselectron, G4bool isrestricted)
83 : fIsRestrictedSamplingRequired(isrestricted), 83 : fIsRestrictedSamplingRequired(isrestricted), fIsElectron(iselectron) {
84 fDCS.resize(gMaxZ+1, nullptr); 84 fDCS.resize(gMaxZ+1, nullptr);
85 fDCSLow.resize(gMaxZ+1, nullptr); 85 fDCSLow.resize(gMaxZ+1, nullptr);
86 fSamplingTables.resize(gMaxZ+1, nullptr); 86 fSamplingTables.resize(gMaxZ+1, nullptr);
87 } 87 }
88 88
89 89
90 // DTR 90 // DTR
91 G4eDPWAElasticDCS::~G4eDPWAElasticDCS() { 91 G4eDPWAElasticDCS::~G4eDPWAElasticDCS() {
92 for (std::size_t i=0; i<fDCS.size(); ++i) { 92 for (std::size_t i=0; i<fDCS.size(); ++i) {
93 if (fDCS[i]) delete fDCS[i]; 93 if (fDCS[i]) delete fDCS[i];
94 } 94 }
95 for (std::size_t i=0; i<fDCSLow.size(); ++i) 95 for (std::size_t i=0; i<fDCSLow.size(); ++i) {
96 if (fDCSLow[i]) delete fDCSLow[i]; 96 if (fDCSLow[i]) delete fDCSLow[i];
97 } 97 }
98 for (std::size_t i=0; i<fSamplingTables.size 98 for (std::size_t i=0; i<fSamplingTables.size(); ++i) {
99 if (fSamplingTables[i]) delete fSamplingTa 99 if (fSamplingTables[i]) delete fSamplingTables[i];
100 } 100 }
101 // clear scp correction data 101 // clear scp correction data
102 for (std::size_t imc=0; imc<fSCPCPerMatCuts. 102 for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
103 if (fSCPCPerMatCuts[imc]) { 103 if (fSCPCPerMatCuts[imc]) {
104 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 104 fSCPCPerMatCuts[imc]->fVSCPC.clear();
105 delete fSCPCPerMatCuts[imc]; 105 delete fSCPCPerMatCuts[imc];
106 } 106 }
107 } 107 }
108 fSCPCPerMatCuts.clear(); 108 fSCPCPerMatCuts.clear();
109 } 109 }
110 110
111 111
112 // initialise for a given 'iz' atomic number: 112 // initialise for a given 'iz' atomic number:
113 // - nothing happens if it has already been i 113 // - nothing happens if it has already been initialised for that Z.
114 void G4eDPWAElasticDCS::InitialiseForZ(std::si 114 void G4eDPWAElasticDCS::InitialiseForZ(std::size_t iz) {
115 if (!gIsGridLoaded) { 115 if (!gIsGridLoaded) {
116 LoadGrid(); 116 LoadGrid();
117 } 117 }
118 LoadDCSForZ((G4int)iz); 118 LoadDCSForZ((G4int)iz);
119 BuildSmplingTableForZ((G4int)iz); 119 BuildSmplingTableForZ((G4int)iz);
120 } 120 }
121 121
122 122
123 // loads the kinetic energy and theta grids fo 123 // loads the kinetic energy and theta grids for the DCS data (first init step)
124 // should be called only by the master 124 // should be called only by the master
125 void G4eDPWAElasticDCS::LoadGrid() { 125 void G4eDPWAElasticDCS::LoadGrid() {
126 G4String fname = FindDirectoryPath() + "grid 126 G4String fname = FindDirectoryPath() + "grid.dat";
127 std::ifstream infile(fname.c_str()); 127 std::ifstream infile(fname.c_str());
128 if (!infile.is_open()) { 128 if (!infile.is_open()) {
129 G4String msg = 129 G4String msg =
130 " Problem while trying to read " + 130 " Problem while trying to read " + fname + " file.\n"+
131 " G4LEDATA version should be G4EML 131 " G4LEDATA version should be G4EMLOW7.12 or later.\n";
132 G4Exception("G4eDPWAElasticDCS::ReadCompre 132 G4Exception("G4eDPWAElasticDCS::ReadCompressedFile","em0006",
133 FatalException,msg.c_str()); 133 FatalException,msg.c_str());
134 return; 134 return;
135 } 135 }
136 // read size 136 // read size
137 infile >> gNumEnergies; 137 infile >> gNumEnergies;
138 infile >> gNumThetas1; 138 infile >> gNumThetas1;
139 infile >> gNumThetas2; 139 infile >> gNumThetas2;
140 // read the grids 140 // read the grids
141 // - energy in [MeV] 141 // - energy in [MeV]
142 G4double dum = 0.0; 142 G4double dum = 0.0;
143 gTheEnergies.resize(gNumEnergies); 143 gTheEnergies.resize(gNumEnergies);
144 for (std::size_t ie=0; ie<gNumEnergies; ++ie 144 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
145 infile >> dum; 145 infile >> dum;
146 gTheEnergies[ie] = G4Log(dum*CLHEP::MeV); 146 gTheEnergies[ie] = G4Log(dum*CLHEP::MeV);
147 if (gTheEnergies[ie]<G4Log(2.0*CLHEP::keV) 147 if (gTheEnergies[ie]<G4Log(2.0*CLHEP::keV)) gIndxEnergyLim = ie; // only for e-
148 } 148 }
149 ++gIndxEnergyLim; 149 ++gIndxEnergyLim;
150 // store/set usefull logarithms of the kinet 150 // store/set usefull logarithms of the kinetic energy grid
151 gLogMinEkin = gTheEnergies[0]; 151 gLogMinEkin = gTheEnergies[0];
152 gInvDelLogEkin = (gNumEnergies-1)/(gTheEnerg 152 gInvDelLogEkin = (gNumEnergies-1)/(gTheEnergies[gNumEnergies-1]-gTheEnergies[0]);
153 // - theta1 in [deg.] (247): we store mu(the 153 // - theta1 in [deg.] (247): we store mu(theta) = 0.5[1-cos(theta)]
154 gTheMus1.resize(gNumThetas1); 154 gTheMus1.resize(gNumThetas1);
155 gTheU1.resize(gNumThetas1); 155 gTheU1.resize(gNumThetas1);
156 const double theA = 0.01; 156 const double theA = 0.01;
157 for (std::size_t it=0; it<gNumThetas1; ++it) 157 for (std::size_t it=0; it<gNumThetas1; ++it) {
158 infile >> dum; 158 infile >> dum;
159 gTheMus1[it] = 0.5*(1.0-std::cos(dum*CLHEP 159 gTheMus1[it] = 0.5*(1.0-std::cos(dum*CLHEP::degree));
160 gTheU1[it] = (theA+1.0)*gTheMus1[it]/(th 160 gTheU1[it] = (theA+1.0)*gTheMus1[it]/(theA+gTheMus1[it]);
161 } 161 }
162 // - theta2 in [deg.] (128): we store mu(the 162 // - theta2 in [deg.] (128): we store mu(theta) = 0.5[1-cos(theta)]
163 gTheMus2.resize(gNumThetas2); 163 gTheMus2.resize(gNumThetas2);
164 gTheU2.resize(gNumThetas2); 164 gTheU2.resize(gNumThetas2);
165 for (std::size_t it=0; it<gNumThetas2; ++it) 165 for (std::size_t it=0; it<gNumThetas2; ++it) {
166 infile >> dum; 166 infile >> dum;
167 gTheMus2[it] = 0.5*(1.0-std::cos(dum*CLHEP 167 gTheMus2[it] = 0.5*(1.0-std::cos(dum*CLHEP::degree));
168 gTheU2[it] = (theA+1.0)*gTheMus2[it]/(th 168 gTheU2[it] = (theA+1.0)*gTheMus2[it]/(theA+gTheMus2[it]);
169 169
170 } 170 }
171 infile.close(); 171 infile.close();
172 gIsGridLoaded = true; 172 gIsGridLoaded = true;
173 } 173 }
174 174
175 175
176 // load DCS data for a given Z 176 // load DCS data for a given Z
177 void G4eDPWAElasticDCS::LoadDCSForZ(G4int iz) 177 void G4eDPWAElasticDCS::LoadDCSForZ(G4int iz) {
178 // Check if it has already been done: 178 // Check if it has already been done:
179 if (fDCS[iz]) return; 179 if (fDCS[iz]) return;
180 // Do it otherwise 180 // Do it otherwise
181 if (fIsElectron) { 181 if (fIsElectron) {
182 // e- 182 // e-
183 // load the high energy part firt: 183 // load the high energy part firt:
184 // - with gNumThetas2 theta and gNumEnergi 184 // - with gNumThetas2 theta and gNumEnergies-gIndxEnergyLim energy values
185 const std::size_t hNumEnergries = gNumEner 185 const std::size_t hNumEnergries = gNumEnergies-gIndxEnergyLim;
186 auto v2DHigh = new G4Physics2DVector(gNumT 186 auto v2DHigh = new G4Physics2DVector(gNumThetas2, hNumEnergries);
187 v2DHigh->SetBicubicInterpolation(true); 187 v2DHigh->SetBicubicInterpolation(true);
188 for (std::size_t it=0; it<gNumThetas2; ++i 188 for (std::size_t it=0; it<gNumThetas2; ++it) {
189 v2DHigh->PutX(it, gTheMus2[it]); 189 v2DHigh->PutX(it, gTheMus2[it]);
190 } 190 }
191 for (std::size_t ie=0; ie<hNumEnergries; + 191 for (std::size_t ie=0; ie<hNumEnergries; ++ie) {
192 v2DHigh->PutY(ie, gTheEnergies[gIndxEner 192 v2DHigh->PutY(ie, gTheEnergies[gIndxEnergyLim+ie]);
193 } 193 }
194 std::ostringstream ossh; 194 std::ostringstream ossh;
195 ossh << FindDirectoryPath() << "dcss/el/dc 195 ossh << FindDirectoryPath() << "dcss/el/dcs_"<< iz<<"_h";
196 std::istringstream finh(std::ios::in); 196 std::istringstream finh(std::ios::in);
197 ReadCompressedFile(ossh.str(), finh); 197 ReadCompressedFile(ossh.str(), finh);
198 G4double dum = 0.0; 198 G4double dum = 0.0;
199 for (std::size_t it=0; it<gNumThetas2; ++i 199 for (std::size_t it=0; it<gNumThetas2; ++it) {
200 finh >> dum; 200 finh >> dum;
201 for (std::size_t ie=0; ie<hNumEnergries; 201 for (std::size_t ie=0; ie<hNumEnergries; ++ie) {
202 finh >> dum; 202 finh >> dum;
203 v2DHigh->PutValue(it, ie, G4Log(dum*CL 203 v2DHigh->PutValue(it, ie, G4Log(dum*CLHEP::cm2/CLHEP::sr));
204 } 204 }
205 } 205 }
206 // load the low energy part: 206 // load the low energy part:
207 // - with gNumThetas1 theta and gIndxEnerg 207 // - with gNumThetas1 theta and gIndxEnergyLim+1 energy values (the +1 is
208 // for including the firts DCS from the 208 // for including the firts DCS from the higher part above for being
209 // able to perform interpolation between 209 // able to perform interpolation between the high and low energy DCS set)
210 auto v2DLow = new G4Physics2DVector(gNumTh 210 auto v2DLow = new G4Physics2DVector(gNumThetas1, gIndxEnergyLim+1);
211 v2DLow->SetBicubicInterpolation(true); 211 v2DLow->SetBicubicInterpolation(true);
212 for (std::size_t it=0; it<gNumThetas1; ++i 212 for (std::size_t it=0; it<gNumThetas1; ++it) {
213 v2DLow->PutX(it, gTheMus1[it]); 213 v2DLow->PutX(it, gTheMus1[it]);
214 } 214 }
215 for (std::size_t ie=0; ie<gIndxEnergyLim+1 215 for (std::size_t ie=0; ie<gIndxEnergyLim+1; ++ie) {
216 v2DLow->PutY(ie, gTheEnergies[ie]); 216 v2DLow->PutY(ie, gTheEnergies[ie]);
217 } 217 }
218 std::ostringstream ossl; 218 std::ostringstream ossl;
219 ossl << FindDirectoryPath() << "dcss/el/dc 219 ossl << FindDirectoryPath() << "dcss/el/dcs_"<< iz<<"_l";
220 std::istringstream finl(std::ios::in); 220 std::istringstream finl(std::ios::in);
221 ReadCompressedFile(ossl.str(), finl); 221 ReadCompressedFile(ossl.str(), finl);
222 for (std::size_t it=0; it<gNumThetas1; ++i 222 for (std::size_t it=0; it<gNumThetas1; ++it) {
223 finl >> dum; 223 finl >> dum;
224 for (std::size_t ie=0; ie<gIndxEnergyLim 224 for (std::size_t ie=0; ie<gIndxEnergyLim; ++ie) {
225 finl >> dum; 225 finl >> dum;
226 v2DLow->PutValue(it, ie, G4Log(dum*CLH 226 v2DLow->PutValue(it, ie, G4Log(dum*CLHEP::cm2/CLHEP::sr));
227 } 227 }
228 } 228 }
229 // add the +1 part: interpolate the firts 229 // add the +1 part: interpolate the firts DCS from the high energy
230 std::size_t ix = 0; 230 std::size_t ix = 0;
231 std::size_t iy = 0; 231 std::size_t iy = 0;
232 for (std::size_t it=0; it<gNumThetas1; ++i 232 for (std::size_t it=0; it<gNumThetas1; ++it) {
233 const G4double val = v2DHigh->Value(gThe 233 const G4double val = v2DHigh->Value(gTheMus1[it], gTheEnergies[gIndxEnergyLim], ix, iy);
234 v2DLow->PutValue(it, gIndxEnergyLim, val 234 v2DLow->PutValue(it, gIndxEnergyLim, val);
235 } 235 }
236 // store 236 // store
237 fDCSLow[iz] = v2DLow; 237 fDCSLow[iz] = v2DLow;
238 fDCS[iz] = v2DHigh; 238 fDCS[iz] = v2DHigh;
239 } else { 239 } else {
240 // e+ 240 // e+
241 auto v2D= new G4Physics2DVector(gNumThetas 241 auto v2D= new G4Physics2DVector(gNumThetas2, gNumEnergies);
242 v2D->SetBicubicInterpolation(true); 242 v2D->SetBicubicInterpolation(true);
243 for (std::size_t it=0; it<gNumThetas2; ++i 243 for (std::size_t it=0; it<gNumThetas2; ++it) {
244 v2D->PutX(it, gTheMus2[it]); 244 v2D->PutX(it, gTheMus2[it]);
245 } 245 }
246 for (std::size_t ie=0; ie<gNumEnergies; ++ 246 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
247 v2D->PutY(ie, gTheEnergies[ie]); 247 v2D->PutY(ie, gTheEnergies[ie]);
248 } 248 }
249 std::ostringstream oss; 249 std::ostringstream oss;
250 oss << FindDirectoryPath() << "dcss/pos/dc 250 oss << FindDirectoryPath() << "dcss/pos/dcs_"<< iz;
251 std::istringstream fin(std::ios::in); 251 std::istringstream fin(std::ios::in);
252 ReadCompressedFile(oss.str(), fin); 252 ReadCompressedFile(oss.str(), fin);
253 G4double dum = 0.0; 253 G4double dum = 0.0;
254 for (std::size_t it=0; it<gNumThetas2; ++i 254 for (std::size_t it=0; it<gNumThetas2; ++it) {
255 fin >> dum; 255 fin >> dum;
256 for (std::size_t ie=0; ie<gNumEnergies; 256 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
257 fin >> dum; 257 fin >> dum;
258 v2D->PutValue(it, ie, G4Log(dum*CLHEP: 258 v2D->PutValue(it, ie, G4Log(dum*CLHEP::cm2/CLHEP::sr));
259 } 259 }
260 } 260 }
261 fDCS[iz]= v2D; 261 fDCS[iz]= v2D;
262 } 262 }
263 } 263 }
264 264
265 265
266 266
267 // Computes the elastic, first and second cros 267 // Computes the elastic, first and second cross sections for the given kinetic
268 // energy and target atom. 268 // energy and target atom.
269 // Cross sections are zero ff ekin is below/ab 269 // Cross sections are zero ff ekin is below/above the kinetic energy grid
270 void G4eDPWAElasticDCS::ComputeCSPerAtom(G4int 270 void G4eDPWAElasticDCS::ComputeCSPerAtom(G4int iz, G4double ekin, G4double& elcs,
271 G4dou 271 G4double& tr1cs, G4double& tr2cs,
272 G4dou 272 G4double mumin, G4double mumax) {
273 // init all cross section values to zero; 273 // init all cross section values to zero;
274 elcs = 0.0; 274 elcs = 0.0;
275 tr1cs = 0.0; 275 tr1cs = 0.0;
276 tr2cs = 0.0; 276 tr2cs = 0.0;
277 // make sure that mu(theta) = 0.5[1-cos(thet 277 // make sure that mu(theta) = 0.5[1-cos(theta)] limits have appropriate vals
278 mumin = std::max(0.0, std::min(1.0, mumin)); 278 mumin = std::max(0.0, std::min(1.0, mumin));
279 mumax = std::max(0.0, std::min(1.0, mumax)); 279 mumax = std::max(0.0, std::min(1.0, mumax));
280 if (mumin>=mumax) return; 280 if (mumin>=mumax) return;
281 // make sure that kin. energy is within the 281 // make sure that kin. energy is within the available range (10 eV-100MeV)
282 const G4double lekin = std::max(gTheEnergies 282 const G4double lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], G4Log(ekin)));
283 // if the lower, denser in theta, DCS set sh 283 // if the lower, denser in theta, DCS set should be used
284 const G4bool isLowerGrid = (fIsElectron && l 284 const G4bool isLowerGrid = (fIsElectron && lekin<gTheEnergies[gIndxEnergyLim]);
285 const std::vector<G4double>& theMuVector = ( 285 const std::vector<G4double>& theMuVector = (isLowerGrid) ? gTheMus1 : gTheMus2;
286 const G4Physics2DVector* the2DDCS = ( 286 const G4Physics2DVector* the2DDCS = (isLowerGrid) ? fDCSLow[iz] : fDCS[iz];
287 // find lower/upper mu bin of integration: 287 // find lower/upper mu bin of integration:
288 // 0.0 <= mumin < 1.0 for sure here 288 // 0.0 <= mumin < 1.0 for sure here
289 const std::size_t iMuStart = (mumin == 0.0) 289 const std::size_t iMuStart = (mumin == 0.0) ? 0 : std::distance( theMuVector.begin(), std::upper_bound(theMuVector.begin(), theMuVector.end(), mumin) )-1 ;
290 // 0.0 < mumax <= 1.0 for sure here 290 // 0.0 < mumax <= 1.0 for sure here
291 const std::size_t iMuEnd = (mumax == 1.0) 291 const std::size_t iMuEnd = (mumax == 1.0) ? theMuVector.size()-2 : std::distance( theMuVector.begin(), std::upper_bound(theMuVector.begin(), theMuVector.end(), mumax) )-1 ;
292 // perform numerical integration of the DCS 292 // perform numerical integration of the DCS over the given [mumin, mumax]
293 // interval (where mu(theta) = 0.5[1-cos(the 293 // interval (where mu(theta) = 0.5[1-cos(theta)]) to get the elastic, first
294 std::size_t ix = 0; 294 std::size_t ix = 0;
295 std::size_t iy = 0; 295 std::size_t iy = 0;
296 for (std::size_t imu=iMuStart; imu<=iMuEnd; 296 for (std::size_t imu=iMuStart; imu<=iMuEnd; ++imu) {
297 G4double elcsPar = 0.0; 297 G4double elcsPar = 0.0;
298 G4double tr1csPar = 0.0; 298 G4double tr1csPar = 0.0;
299 G4double tr2csPar = 0.0; 299 G4double tr2csPar = 0.0;
300 const G4double low = (imu==iMuStart) ? mum 300 const G4double low = (imu==iMuStart) ? mumin : theMuVector[imu];
301 const G4double del = (imu==iMuEnd) ? mum 301 const G4double del = (imu==iMuEnd) ? mumax-low : theMuVector[imu+1]-low;
302 ix = imu; 302 ix = imu;
303 for (std::size_t igl=0; igl<8; ++igl) { 303 for (std::size_t igl=0; igl<8; ++igl) {
304 const double mu = low + del*gXGL[igl]; 304 const double mu = low + del*gXGL[igl];
305 const double dcs = G4Exp(the2DDCS->Value 305 const double dcs = G4Exp(the2DDCS->Value(mu, lekin, ix, iy));
306 elcsPar += gWGL[igl]*dcs; / 306 elcsPar += gWGL[igl]*dcs; // elastic
307 tr1csPar += gWGL[igl]*dcs*mu; / 307 tr1csPar += gWGL[igl]*dcs*mu; // first transport
308 tr2csPar += gWGL[igl]*dcs*mu*(1.0-mu); / 308 tr2csPar += gWGL[igl]*dcs*mu*(1.0-mu); // second transport
309 } 309 }
310 elcs += del*elcsPar; 310 elcs += del*elcsPar;
311 tr1cs += del*tr1csPar; 311 tr1cs += del*tr1csPar;
312 tr2cs += del*tr2csPar; 312 tr2cs += del*tr2csPar;
313 } 313 }
314 elcs *= 2.0*CLHEP::twopi; 314 elcs *= 2.0*CLHEP::twopi;
315 tr1cs *= 4.0*CLHEP::twopi; 315 tr1cs *= 4.0*CLHEP::twopi;
316 tr2cs *= 12.0*CLHEP::twopi; 316 tr2cs *= 12.0*CLHEP::twopi;
317 } 317 }
318 318
319 319
320 // data structure to store one sampling table: 320 // data structure to store one sampling table: combined Alias + RatIn
321 // NOTE: when Alias is used, sampling on a res 321 // NOTE: when Alias is used, sampling on a resctricted interval is not possible
322 // However, Alias makes possible faster 322 // However, Alias makes possible faster sampling. Alias is used in case
323 // of single scattering model while it's 323 // of single scattering model while it's not used in case of mixed-model
324 // when restricted interval sampling is 324 // when restricted interval sampling is needed. This is controlled by
325 // the fIsRestrictedSamplingRequired fla 325 // the fIsRestrictedSamplingRequired flag (false by default).
326 struct OneSamplingTable { 326 struct OneSamplingTable {
327 OneSamplingTable () = default; 327 OneSamplingTable () = default;
328 void SetSize(std::size_t nx, G4bool useAlias 328 void SetSize(std::size_t nx, G4bool useAlias) {
329 fN = nx; 329 fN = nx;
330 // Alias 330 // Alias
331 if (useAlias) { 331 if (useAlias) {
332 fW.resize(nx); 332 fW.resize(nx);
333 fI.resize(nx); 333 fI.resize(nx);
334 } 334 }
335 // Ratin 335 // Ratin
336 fCum.resize(nx); 336 fCum.resize(nx);
337 fA.resize(nx); 337 fA.resize(nx);
338 fB.resize(nx); 338 fB.resize(nx);
339 } 339 }
340 340
341 // members 341 // members
342 std::size_t fN; // # da 342 std::size_t fN; // # data points
343 G4double fScreenParA; // the 343 G4double fScreenParA; // the screening parameter
344 std::vector<G4double> fW; 344 std::vector<G4double> fW;
345 std::vector<G4double> fCum; 345 std::vector<G4double> fCum;
346 std::vector<G4double> fA; 346 std::vector<G4double> fA;
347 std::vector<G4double> fB; 347 std::vector<G4double> fB;
348 std::vector<G4int> fI; 348 std::vector<G4int> fI;
349 }; 349 };
350 350
351 351
352 // loads sampling table for the given Z over t 352 // loads sampling table for the given Z over the enrgy grid
353 void G4eDPWAElasticDCS::BuildSmplingTableForZ( 353 void G4eDPWAElasticDCS::BuildSmplingTableForZ(G4int iz) {
354 // Check if it has already been done: 354 // Check if it has already been done:
355 if (fSamplingTables[iz]) return; 355 if (fSamplingTables[iz]) return;
356 // Do it otherwise: 356 // Do it otherwise:
357 // allocate space 357 // allocate space
358 auto sTables = new std::vector<OneSamplingTa 358 auto sTables = new std::vector<OneSamplingTable>(gNumEnergies);
359 // read compressed sampling table data 359 // read compressed sampling table data
360 std::ostringstream oss; 360 std::ostringstream oss;
361 const G4String fname = fIsElectron ? "stable 361 const G4String fname = fIsElectron ? "stables/el/" : "stables/pos/";
362 oss << FindDirectoryPath() << fname << "stab 362 oss << FindDirectoryPath() << fname << "stable_" << iz;
363 std::istringstream fin(std::ios::in); 363 std::istringstream fin(std::ios::in);
364 ReadCompressedFile(oss.str(), fin); 364 ReadCompressedFile(oss.str(), fin);
365 std::size_t ndata = 0; 365 std::size_t ndata = 0;
366 for (std::size_t ie=0; ie<gNumEnergies; ++ie 366 for (std::size_t ie=0; ie<gNumEnergies; ++ie) {
367 OneSamplingTable& aTable = (*sTables)[ie]; 367 OneSamplingTable& aTable = (*sTables)[ie];
368 // #data in this table 368 // #data in this table
369 fin >> ndata; 369 fin >> ndata;
370 aTable.SetSize(ndata, !fIsRestrictedSampli 370 aTable.SetSize(ndata, !fIsRestrictedSamplingRequired);
371 // the A screening parameter value used fo 371 // the A screening parameter value used for transformation of mu to u
372 fin >> aTable.fScreenParA; 372 fin >> aTable.fScreenParA;
373 // load data: Alias(W,I) + RatIn(Cum, A, B 373 // load data: Alias(W,I) + RatIn(Cum, A, B)
374 if (!fIsRestrictedSamplingRequired) { 374 if (!fIsRestrictedSamplingRequired) {
375 for (std::size_t id=0; id<ndata; ++id) { 375 for (std::size_t id=0; id<ndata; ++id) {
376 fin >> aTable.fW[id]; 376 fin >> aTable.fW[id];
377 } 377 }
378 for (std::size_t id=0; id<ndata; ++id) { 378 for (std::size_t id=0; id<ndata; ++id) {
379 fin >> aTable.fI[id]; 379 fin >> aTable.fI[id];
380 } 380 }
381 } 381 }
382 for (std::size_t id=0; id<ndata; ++id) { 382 for (std::size_t id=0; id<ndata; ++id) {
383 fin >> aTable.fCum[id]; 383 fin >> aTable.fCum[id];
384 } 384 }
385 for (std::size_t id=0; id<ndata; ++id) { 385 for (std::size_t id=0; id<ndata; ++id) {
386 fin >> aTable.fA[id]; 386 fin >> aTable.fA[id];
387 } 387 }
388 for (std::size_t id=0; id<ndata; ++id) { 388 for (std::size_t id=0; id<ndata; ++id) {
389 fin >> aTable.fB[id]; 389 fin >> aTable.fB[id];
390 } 390 }
391 } 391 }
392 fSamplingTables[iz] = sTables; 392 fSamplingTables[iz] = sTables;
393 } 393 }
394 394
395 395
396 396
397 397
398 398
399 // samples cos(theta) i.e. cosine of the polar 399 // samples cos(theta) i.e. cosine of the polar angle of scattering in elastic
400 // interaction (Coulomb scattering) of the pro 400 // interaction (Coulomb scattering) of the projectile (e- or e+ depending on
401 // fIsElectron) with kinetic energy of exp('le 401 // fIsElectron) with kinetic energy of exp('lekin'), target atom with atomic
402 // muber of 'iz'. See the 'SampleCosineThetaRe 402 // muber of 'iz'. See the 'SampleCosineThetaRestricted' for obtain samples on
403 // a restricted inteval. 403 // a restricted inteval.
404 G4double 404 G4double
405 G4eDPWAElasticDCS::SampleCosineTheta(std::size 405 G4eDPWAElasticDCS::SampleCosineTheta(std::size_t iz, G4double lekin, G4double r1,
406 G4double 406 G4double r2, G4double r3) {
407 lekin = std::max(gTheEnergies[0], std::min(g 407 lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], lekin));
408 // determine the discrete ekin sampling tabl 408 // determine the discrete ekin sampling table to be used:
409 // - statistical interpolation (i.e. linear 409 // - statistical interpolation (i.e. linear) on log energy scale
410 const G4double rem = (lekin-gLogMinEkin 410 const G4double rem = (lekin-gLogMinEkin)*gInvDelLogEkin;
411 const auto k = (std::size_t)rem; 411 const auto k = (std::size_t)rem;
412 const std::size_t iekin = (r1 < rem-k) ? k+1 412 const std::size_t iekin = (r1 < rem-k) ? k+1 : k;
413 // sample the mu(t)=0.5(1-cos(t)) 413 // sample the mu(t)=0.5(1-cos(t))
414 const double mu = SampleMu(iz, iekin, r2, 414 const double mu = SampleMu(iz, iekin, r2, r3);
415 return std::max(-1.0, std::min(1.0, 1.0-2.0* 415 return std::max(-1.0, std::min(1.0, 1.0-2.0*mu));
416 } 416 }
417 417
418 418
419 // samples cos(theta) i.e. cosine of the polar 419 // samples cos(theta) i.e. cosine of the polar angle of scattering in elastic
420 // interaction (Coulomb scattering) of the pro 420 // interaction (Coulomb scattering) of the projectile (e- or e+ depending on
421 // fIsElectron) with kinetic energy of exp('le 421 // fIsElectron) with kinetic energy of exp('lekin'), target atom with atomic
422 // muber of 'iz'. 422 // muber of 'iz'.
423 // The cosine theta will be in the [costMin, c 423 // The cosine theta will be in the [costMin, costMax] interval where costMin
424 // corresponds to a maximum allowed polar scat 424 // corresponds to a maximum allowed polar scattering angle thetaMax while
425 // costMin corresponds to minimum allowed pola 425 // costMin corresponds to minimum allowed polar scatterin angle thetaMin.
426 // See the 'SampleCosineTheta' for obtain samp 426 // See the 'SampleCosineTheta' for obtain samples on the entire [-1,1] range.
427 G4double 427 G4double
428 G4eDPWAElasticDCS::SampleCosineThetaRestricted 428 G4eDPWAElasticDCS::SampleCosineThetaRestricted(std::size_t iz, G4double lekin,
429 429 G4double r1, G4double r2,
430 430 G4double costMax, G4double costMin) {
431 // costMin corresponds to mu-max while costM 431 // costMin corresponds to mu-max while costMax to mu-min: mu(t)=0.5[1-cos(t)]
432 lekin = std::max(gTheEnergies[0], std::min(g 432 lekin = std::max(gTheEnergies[0], std::min(gTheEnergies[gNumEnergies-1], lekin));
433 // determine the discrete ekin sampling tabl 433 // determine the discrete ekin sampling table to be used:
434 // - statistical interpolation (i.e. linear 434 // - statistical interpolation (i.e. linear) on log energy scale
435 const G4double rem = (lekin-gLogMinEkin 435 const G4double rem = (lekin-gLogMinEkin)*gInvDelLogEkin;
436 const auto k = (size_t)rem; 436 const auto k = (size_t)rem;
437 const std::size_t iekin = (r1 < rem-k) ? k : 437 const std::size_t iekin = (r1 < rem-k) ? k : k+1;
438 // sample the mu(t)=0.5(1-cos(t)) 438 // sample the mu(t)=0.5(1-cos(t))
439 const G4double mu = SampleMu(iz, iekin, r2, 439 const G4double mu = SampleMu(iz, iekin, r2, 0.5*(1.0-costMax), 0.5*(1.0-costMin));
440 return std::max(-1.0, std::min(1.0, 1.0-2.0* 440 return std::max(-1.0, std::min(1.0, 1.0-2.0*mu));
441 } 441 }
442 442
443 443
444 444
445 G4double 445 G4double
446 G4eDPWAElasticDCS::SampleMu(std::size_t izet, 446 G4eDPWAElasticDCS::SampleMu(std::size_t izet, std::size_t ie, G4double r1, G4double r2) {
447 OneSamplingTable& rtn = (*fSamplingTables[iz 447 OneSamplingTable& rtn = (*fSamplingTables[izet])[ie];
448 // get the lower index of the bin by using t 448 // get the lower index of the bin by using the alias part
449 const G4double rest = r1 * (rtn.fN - 1); 449 const G4double rest = r1 * (rtn.fN - 1);
450 auto indxl = (std::size_t)rest; 450 auto indxl = (std::size_t)rest;
451 const G4double dum0 = rest - indxl; 451 const G4double dum0 = rest - indxl;
452 if (rtn.fW[indxl] < dum0) indxl = rtn.fI[ind 452 if (rtn.fW[indxl] < dum0) indxl = rtn.fI[indxl];
453 // sample value within the selected bin by u 453 // sample value within the selected bin by using ratin based numerical inversion
454 const G4double delta = rtn.fCum[indxl + 1] - 454 const G4double delta = rtn.fCum[indxl + 1] - rtn.fCum[indxl];
455 const G4double aval = r2 * delta; 455 const G4double aval = r2 * delta;
456 456
457 const G4double dum1 = (1.0 + rtn.fA[indxl] 457 const G4double dum1 = (1.0 + rtn.fA[indxl] + rtn.fB[indxl]) * delta * aval;
458 const G4double dum2 = delta * delta + rtn.f 458 const G4double dum2 = delta * delta + rtn.fA[indxl] * delta * aval + rtn.fB[indxl] * aval * aval;
459 const std::vector<G4double>& theUVect = (fIs 459 const std::vector<G4double>& theUVect = (fIsElectron && ie<gIndxEnergyLim) ? gTheU1 : gTheU2;
460 const G4double u = theUVect[indxl] + dum 460 const G4double u = theUVect[indxl] + dum1 / dum2 * (theUVect[indxl + 1] - theUVect[indxl]);
461 // transform back u to mu 461 // transform back u to mu
462 return rtn.fScreenParA*u/(rtn.fScreenParA+1. 462 return rtn.fScreenParA*u/(rtn.fScreenParA+1.0-u);
463 } 463 }
464 464
465 465
466 466
467 G4double 467 G4double
468 G4eDPWAElasticDCS::FindCumValue(G4double u, co 468 G4eDPWAElasticDCS::FindCumValue(G4double u, const OneSamplingTable& stable,
469 const std::vec 469 const std::vector<G4double>& uvect) {
470 const std::size_t iLow = std::distance( uvec 470 const std::size_t iLow = std::distance( uvect.begin(), std::upper_bound(uvect.begin(), uvect.end(), u) )-1;
471 const G4double tau = (u-uvect[iLow])/(uv 471 const G4double tau = (u-uvect[iLow])/(uvect[iLow+1]-uvect[iLow]); // Note: I could store 1/(fX[iLow+1]-fX[iLow])
472 const G4double dum0 = (1.0+stable.fA[iLow 472 const G4double dum0 = (1.0+stable.fA[iLow]*(1.0-tau)+stable.fB[iLow]);
473 const G4double dum1 = 2.0*stable.fB[iLow] 473 const G4double dum1 = 2.0*stable.fB[iLow]*tau;
474 const G4double dum2 = 1.0 - std::sqrt(std 474 const G4double dum2 = 1.0 - std::sqrt(std::max(0.0, 1.0-2.0*dum1*tau/(dum0*dum0)));
475 return std::min(stable.fCum[iLow+1], std::ma 475 return std::min(stable.fCum[iLow+1], std::max(stable.fCum[iLow], stable.fCum[iLow]+dum0*dum2*(stable.fCum[iLow+1]-stable.fCum[iLow])/dum1 ));
476 } 476 }
477 477
478 // muMin and muMax : no checks on these 478 // muMin and muMax : no checks on these
479 G4double G4eDPWAElasticDCS::SampleMu(std::size 479 G4double G4eDPWAElasticDCS::SampleMu(std::size_t izet, std::size_t ie, G4double r1,
480 G4double 480 G4double muMin, G4double muMax) {
481 const OneSamplingTable& rtn = (*fSamplingTab 481 const OneSamplingTable& rtn = (*fSamplingTables[izet])[ie];
482 const G4double theA = rtn.fScreenParA; 482 const G4double theA = rtn.fScreenParA;
483 // 483 //
484 const std::vector<G4double>& theUVect = (fIs 484 const std::vector<G4double>& theUVect = (fIsElectron && ie<gIndxEnergyLim) ? gTheU1 : gTheU2;
485 const G4double xiMin = (muMin > 0.0) ? Fin 485 const G4double xiMin = (muMin > 0.0) ? FindCumValue((theA+1.0)*muMin/(theA+muMin), rtn, theUVect) : 0.0;
486 const G4double xiMax = (muMax < 1.0) ? Fin 486 const G4double xiMax = (muMax < 1.0) ? FindCumValue((theA+1.0)*muMax/(theA+muMax), rtn, theUVect) : 1.0;
487 // 487 //
488 const G4double xi = xiMin+r1*(xiMax-xiM 488 const G4double xi = xiMin+r1*(xiMax-xiMin); // a smaple within the range
489 const std::size_t iLow = std::distance( rtn. 489 const std::size_t iLow = std::distance( rtn.fCum.begin(), std::upper_bound(rtn.fCum.begin(), rtn.fCum.end(), xi) )-1;
490 const G4double delta = rtn.fCum[iLow + 1] 490 const G4double delta = rtn.fCum[iLow + 1] - rtn.fCum[iLow];
491 const G4double aval = xi - rtn.fCum[iLow] 491 const G4double aval = xi - rtn.fCum[iLow];
492 492
493 const G4double dum1 = (1.0 + rtn.fA[iLow] 493 const G4double dum1 = (1.0 + rtn.fA[iLow] + rtn.fB[iLow]) * delta * aval;
494 const G4double dum2 = delta * delta + rtn 494 const G4double dum2 = delta * delta + rtn.fA[iLow] * delta * aval + rtn.fB[iLow] * aval * aval;
495 const G4double u = theUVect[iLow] + du 495 const G4double u = theUVect[iLow] + dum1 / dum2 * (theUVect[iLow + 1] - theUVect[iLow]);
496 return theA*u/(theA+1.0-u); 496 return theA*u/(theA+1.0-u);
497 } 497 }
498 498
499 499
500 500
501 // set the DCS data directory path 501 // set the DCS data directory path
502 const G4String& G4eDPWAElasticDCS::FindDirecto 502 const G4String& G4eDPWAElasticDCS::FindDirectoryPath() {
503 // check environment variable 503 // check environment variable
504 if (gDataDirectory.empty()) { 504 if (gDataDirectory.empty()) {
505 std::ostringstream ost; 505 std::ostringstream ost;
506 ost << G4EmParameters::Instance()->GetDirL 506 ost << G4EmParameters::Instance()->GetDirLEDATA() << "/dpwa/";
507 gDataDirectory = ost.str(); 507 gDataDirectory = ost.str();
508 } 508 }
509 return gDataDirectory; 509 return gDataDirectory;
510 } 510 }
511 511
512 512
513 513
514 // uncompress one data file into the input str 514 // uncompress one data file into the input string stream
515 void 515 void
516 G4eDPWAElasticDCS::ReadCompressedFile(G4String 516 G4eDPWAElasticDCS::ReadCompressedFile(G4String fname, std::istringstream &iss) {
517 G4String *dataString = nullptr; 517 G4String *dataString = nullptr;
518 G4String compfilename(fname+".z"); 518 G4String compfilename(fname+".z");
519 // create input stream with binary mode oper 519 // create input stream with binary mode operation and positioning at the end of the file
520 std::ifstream in(compfilename, std::ios::bin 520 std::ifstream in(compfilename, std::ios::binary | std::ios::ate);
521 if (in.good()) { 521 if (in.good()) {
522 // get current position in the stream (wa 522 // get current position in the stream (was set to the end)
523 std::size_t fileSize = in.tellg(); 523 std::size_t fileSize = in.tellg();
524 // set current position being the beginni 524 // set current position being the beginning of the stream
525 in.seekg(0,std::ios::beg); 525 in.seekg(0,std::ios::beg);
526 // create (zlib) byte buffer for the data 526 // create (zlib) byte buffer for the data
527 Bytef *compdata = new Bytef[fileSize]; 527 Bytef *compdata = new Bytef[fileSize];
528 while(in) { 528 while(in) {
529 in.read((char*)compdata, fileSize); 529 in.read((char*)compdata, fileSize);
530 } 530 }
531 // create (zlib) byte buffer for the unco 531 // create (zlib) byte buffer for the uncompressed data
532 uLongf complen = (uLongf)(fileSize*4); 532 uLongf complen = (uLongf)(fileSize*4);
533 Bytef *uncompdata = new Bytef[complen]; 533 Bytef *uncompdata = new Bytef[complen];
534 while (Z_OK!=uncompress(uncompdata, &comp 534 while (Z_OK!=uncompress(uncompdata, &complen, compdata, fileSize)) {
535 // increase uncompressed byte buffer 535 // increase uncompressed byte buffer
536 delete[] uncompdata; 536 delete[] uncompdata;
537 complen *= 2; 537 complen *= 2;
538 uncompdata = new Bytef[complen]; 538 uncompdata = new Bytef[complen];
539 } 539 }
540 // delete the compressed data buffer 540 // delete the compressed data buffer
541 delete [] compdata; 541 delete [] compdata;
542 // create a string from the uncompressed 542 // create a string from the uncompressed data (will be deallocated by the caller)
543 dataString = new G4String((char*)uncompda 543 dataString = new G4String((char*)uncompdata, (long)complen);
544 // delete the uncompressed data buffer 544 // delete the uncompressed data buffer
545 delete [] uncompdata; 545 delete [] uncompdata;
546 } else { 546 } else {
547 G4String msg = 547 G4String msg =
548 " Problem while trying to read " + 548 " Problem while trying to read " + fname + " data file.\n"+
549 " G4LEDATA version should be G4EML 549 " G4LEDATA version should be G4EMLOW7.12 or later.\n";
550 G4Exception("G4eDPWAElasticDCS::ReadCompre 550 G4Exception("G4eDPWAElasticDCS::ReadCompressedFile","em0006",
551 FatalException,msg.c_str()); 551 FatalException,msg.c_str());
552 return; 552 return;
553 } 553 }
554 // create the input string stream from the d 554 // create the input string stream from the data string
555 if (dataString) { 555 if (dataString) {
556 iss.str(*dataString); 556 iss.str(*dataString);
557 in.close(); 557 in.close();
558 delete dataString; 558 delete dataString;
559 } 559 }
560 } 560 }
561 561
562 562
563 563
564 564
565 G4double 565 G4double
566 G4eDPWAElasticDCS::ComputeScatteringPowerCorre 566 G4eDPWAElasticDCS::ComputeScatteringPowerCorrection(const G4MaterialCutsCouple *matcut,
567 567 G4double ekin) {
568 const G4int imc = matcut->GetIndex(); 568 const G4int imc = matcut->GetIndex();
569 G4double corFactor = 1.0; 569 G4double corFactor = 1.0;
570 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin< 570 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<=fSCPCPerMatCuts[imc]->fPrCut) {
571 return corFactor; 571 return corFactor;
572 } 572 }
573 // get the scattering power correction facto 573 // get the scattering power correction factor
574 const G4double lekin = G4Log(ekin); 574 const G4double lekin = G4Log(ekin);
575 G4double remaining = (lekin-fSCPCPerMatCut 575 G4double remaining = (lekin-fSCPCPerMatCuts[imc]->fLEmin)*fSCPCPerMatCuts[imc]->fILDel;
576 std::size_t lindx = (G4int)remaining; 576 std::size_t lindx = (G4int)remaining;
577 remaining -= lindx; 577 remaining -= lindx;
578 std::size_t imax = fSCPCPerMatCuts[imc]- 578 std::size_t imax = fSCPCPerMatCuts[imc]->fVSCPC.size()-1;
579 if (lindx>=imax) { 579 if (lindx>=imax) {
580 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i 580 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[imax];
581 } else { 581 } else {
582 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l 582 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[lindx] + remaining*(fSCPCPerMatCuts[imc]->fVSCPC[lindx+1]-fSCPCPerMatCuts[imc]->fVSCPC[lindx]);
583 } 583 }
584 return corFactor; 584 return corFactor;
585 } 585 }
586 586
587 587
588 void G4eDPWAElasticDCS::InitSCPCorrection(G4do 588 void G4eDPWAElasticDCS::InitSCPCorrection(G4double lowEnergyLimit,
589 G4do 589 G4double highEnergyLimit) {
590 // get the material-cuts table 590 // get the material-cuts table
591 G4ProductionCutsTable *thePCTable = G4Produc 591 G4ProductionCutsTable *thePCTable = G4ProductionCutsTable::GetProductionCutsTable();
592 std::size_t numMatCuts = thePCTab 592 std::size_t numMatCuts = thePCTable->GetTableSize();
593 // clear container if any 593 // clear container if any
594 for (std::size_t imc=0; imc<fSCPCPerMatCuts. 594 for (std::size_t imc=0; imc<fSCPCPerMatCuts.size(); ++imc) {
595 if (fSCPCPerMatCuts[imc]) { 595 if (fSCPCPerMatCuts[imc]) {
596 fSCPCPerMatCuts[imc]->fVSCPC.clear(); 596 fSCPCPerMatCuts[imc]->fVSCPC.clear();
597 delete fSCPCPerMatCuts[imc]; 597 delete fSCPCPerMatCuts[imc];
598 fSCPCPerMatCuts[imc] = nullptr; 598 fSCPCPerMatCuts[imc] = nullptr;
599 } 599 }
600 } 600 }
601 // 601 //
602 // set size of the container and create the 602 // set size of the container and create the corresponding data structures
603 fSCPCPerMatCuts.resize(numMatCuts,nullptr); 603 fSCPCPerMatCuts.resize(numMatCuts,nullptr);
604 // loop over the material-cuts and create sc 604 // loop over the material-cuts and create scattering power correction data structure for each
605 for (G4int imc=0; imc<(G4int)numMatCuts; ++i 605 for (G4int imc=0; imc<(G4int)numMatCuts; ++imc) {
606 const G4MaterialCutsCouple *matCut = theP 606 const G4MaterialCutsCouple *matCut = thePCTable->GetMaterialCutsCouple(imc);
607 const G4Material* mat = matCut->GetMateri 607 const G4Material* mat = matCut->GetMaterial();
608 // get e- production cut in the current ma 608 // get e- production cut in the current material-cuts in energy
609 const G4double ecut = (*(thePCTable->Ge 609 const G4double ecut = (*(thePCTable->GetEnergyCutsVector(idxG4ElectronCut)))[matCut->GetIndex()];
610 const G4double limit = fIsElectron ? 2.0 610 const G4double limit = fIsElectron ? 2.0*ecut : ecut;
611 const G4double min = std::max(limit,lo 611 const G4double min = std::max(limit,lowEnergyLimit);
612 const G4double max = highEnergyLimit; 612 const G4double max = highEnergyLimit;
613 if (min>=max) { 613 if (min>=max) {
614 fSCPCPerMatCuts[imc] = new SCPCorrection 614 fSCPCPerMatCuts[imc] = new SCPCorrection();
615 fSCPCPerMatCuts[imc]->fIsUse = false; 615 fSCPCPerMatCuts[imc]->fIsUse = false;
616 fSCPCPerMatCuts[imc]->fPrCut = min; 616 fSCPCPerMatCuts[imc]->fPrCut = min;
617 continue; 617 continue;
618 } 618 }
619 G4int numEbins = fNumSPCEbinPerDec 619 G4int numEbins = fNumSPCEbinPerDec*G4lrint(std::log10(max/min));
620 numEbins = std::max(numEbins 620 numEbins = std::max(numEbins,3);
621 const G4double lmin = G4Log(min); 621 const G4double lmin = G4Log(min);
622 const G4double ldel = G4Log(max/min)/(n 622 const G4double ldel = G4Log(max/min)/(numEbins-1.0);
623 fSCPCPerMatCuts[imc] = new SCPCorrection 623 fSCPCPerMatCuts[imc] = new SCPCorrection();
624 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi 624 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbins,1.0);
625 fSCPCPerMatCuts[imc]->fIsUse = true; 625 fSCPCPerMatCuts[imc]->fIsUse = true;
626 fSCPCPerMatCuts[imc]->fPrCut = min; 626 fSCPCPerMatCuts[imc]->fPrCut = min;
627 fSCPCPerMatCuts[imc]->fLEmin = lmin; 627 fSCPCPerMatCuts[imc]->fLEmin = lmin;
628 fSCPCPerMatCuts[imc]->fILDel = 1./ldel; 628 fSCPCPerMatCuts[imc]->fILDel = 1./ldel;
629 // compute Moliere material dependet param 629 // compute Moliere material dependet parameetrs
630 G4double moliereBc = 0.0; 630 G4double moliereBc = 0.0;
631 G4double moliereXc2 = 0.0; 631 G4double moliereXc2 = 0.0;
632 ComputeMParams(mat, moliereBc, moliereXc2) 632 ComputeMParams(mat, moliereBc, moliereXc2);
633 // compute scattering power correction ove 633 // compute scattering power correction over the enrgy grid
634 for (G4int ie=0; ie<numEbins; ++ie) { 634 for (G4int ie=0; ie<numEbins; ++ie) {
635 const G4double ekin = G4Exp(lmin+ie*ldel 635 const G4double ekin = G4Exp(lmin+ie*ldel);
636 G4double scpCorr = 1.0; 636 G4double scpCorr = 1.0;
637 // compute correction factor: I.Kawrakow 637 // compute correction factor: I.Kawrakow, Med.Phys.24,505-517(1997)(Eqs(32-37)
638 if (ie>0) { 638 if (ie>0) {
639 const G4double tau = ekin/CLHEP:: 639 const G4double tau = ekin/CLHEP::electron_mass_c2;
640 const G4double tauCut = ecut/CLHEP:: 640 const G4double tauCut = ecut/CLHEP::electron_mass_c2;
641 // Moliere's screening parameter 641 // Moliere's screening parameter
642 const G4double A = moliereXc2/( 642 const G4double A = moliereXc2/(4.0*tau*(tau+2.)*moliereBc);
643 const G4double gr = (1.+2.*A)*G4 643 const G4double gr = (1.+2.*A)*G4Log(1.+1./A)-2.;
644 const G4double dum0 = (tau+2.)/(ta 644 const G4double dum0 = (tau+2.)/(tau+1.);
645 const G4double dum1 = tau+1.; 645 const G4double dum1 = tau+1.;
646 G4double gm = G4Log(0.5*tau/tauCut) 646 G4double gm = G4Log(0.5*tau/tauCut) + (1.+dum0*dum0)*G4Log(2.*(tau-tauCut+2.)/(tau+4.))
647 - 0.25*(tau+2.)*( tau+ 647 - 0.25*(tau+2.)*( tau+2.+2.*(2.*tau+1.)/(dum1*dum1))*
648 G4Log((tau+4.)*(tau- 648 G4Log((tau+4.)*(tau-tauCut)/tau/(tau-tauCut+2.))
649 + 0.5*(tau-2*tauCut)*( 649 + 0.5*(tau-2*tauCut)*(tau+2.)*(1./(tau-tauCut)-1./(dum1*dum1));
650 if (gm<gr) { 650 if (gm<gr) {
651 gm = gm/gr; 651 gm = gm/gr;
652 } else { 652 } else {
653 gm = 1.; 653 gm = 1.;
654 } 654 }
655 const G4double z0 = matCut->GetMateri 655 const G4double z0 = matCut->GetMaterial()->GetIonisation()->GetZeffective();
656 scpCorr = 1.-gm*z0/(z0*(z0+1.)); 656 scpCorr = 1.-gm*z0/(z0*(z0+1.));
657 } 657 }
658 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo 658 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCorr;
659 } 659 }
660 } 660 }
661 } 661 }
662 662
663 663
664 // compute material dependent Moliere MSC para 664 // compute material dependent Moliere MSC parameters at initialisation
665 void G4eDPWAElasticDCS::ComputeMParams(const G 665 void G4eDPWAElasticDCS::ComputeMParams(const G4Material* mat, G4double& theBc,
666 G4doubl 666 G4double& theXc2) {
667 const G4double const1 = 7821.6; // [ 667 const G4double const1 = 7821.6; // [cm2/g]
668 const G4double const2 = 0.1569; // [ 668 const G4double const2 = 0.1569; // [cm2 MeV2 / g]
669 const G4double finstrc2 = 5.325135453E-5; / 669 const G4double finstrc2 = 5.325135453E-5; // fine-structure const. square
670 // G4double xi = 1.0; 670 // G4double xi = 1.0;
671 const G4ElementVector* theElemVect = ma 671 const G4ElementVector* theElemVect = mat->GetElementVector();
672 const std::size_t numelems = ma 672 const std::size_t numelems = mat->GetNumberOfElements();
673 // 673 //
674 const G4double* theNbAtomsPerVolVect = ma 674 const G4double* theNbAtomsPerVolVect = mat->GetVecNbOfAtomsPerVolume();
675 G4double theTotNbAtomsPerVol = ma 675 G4double theTotNbAtomsPerVol = mat->GetTotNbOfAtomsPerVolume();
676 // 676 //
677 G4double zs = 0.0; 677 G4double zs = 0.0;
678 G4double zx = 0.0; 678 G4double zx = 0.0;
679 G4double ze = 0.0; 679 G4double ze = 0.0;
680 G4double sa = 0.0; 680 G4double sa = 0.0;
681 // 681 //
682 for(std::size_t ielem = 0; ielem < numelems 682 for(std::size_t ielem = 0; ielem < numelems; ++ielem) {
683 const G4double zet = (*theElemVect)[iele 683 const G4double zet = (*theElemVect)[ielem]->GetZ();
684 const G4double iwa = (*theElemVect)[iele 684 const G4double iwa = (*theElemVect)[ielem]->GetN();
685 const G4double ipz = theNbAtomsPerVolVec 685 const G4double ipz = theNbAtomsPerVolVect[ielem]/theTotNbAtomsPerVol;
686 const G4double dum = ipz*zet*(zet+1.0); 686 const G4double dum = ipz*zet*(zet+1.0);
687 zs += dum; 687 zs += dum;
688 ze += dum*(-2.0/3.0)*G4Log(zet) 688 ze += dum*(-2.0/3.0)*G4Log(zet);
689 zx += dum*G4Log(1.0+3.34*finstr 689 zx += dum*G4Log(1.0+3.34*finstrc2*zet*zet);
690 sa += ipz*iwa; 690 sa += ipz*iwa;
691 } 691 }
692 const G4double density = mat->GetDensity()* 692 const G4double density = mat->GetDensity()*CLHEP::cm3/CLHEP::g; // [g/cm3]
693 // 693 //
694 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa; << 694 theBc = const1*density*zs/sa*G4Exp(ze/zs)/G4Exp(zx/zs); //[1/cm]
695 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx) << 695 theXc2 = const2*density*zs/sa; // [MeV2/cm]
696 <<
697 theBc = const1*density*z0*G4Exp(z1); //[1 <<
698 theXc2 = const2*density*z0; // [MeV2/cm] <<
699 // change to Geant4 internal units of 1/len 696 // change to Geant4 internal units of 1/length and energ2/length
700 theBc *= 1.0/CLHEP::cm; 697 theBc *= 1.0/CLHEP::cm;
701 theXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm; 698 theXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
702 } 699 }
703 700
704 701