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