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