| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
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::gIsGr
53 G4String G4eDPWAElasticDCS::gData
54 // final values of these variables will be set
55 std::size_t G4eDPWAElasticDCS::gNumE
56 std::size_t G4eDPWAElasticDCS::gIndx
57 std::size_t G4eDPWAElasticDCS::gNumT
58 std::size_t G4eDPWAElasticDCS::gNumT
59 G4double G4eDPWAElasticDCS::gLogM
60 G4double G4eDPWAElasticDCS::gInvD
61 // containers for grids: Ekin, mu(t)=0.5[1-cos
62 std::vector<G4double> G4eDPWAElasticDCS::gTheE
63 std::vector<G4double> G4eDPWAElasticDCS::gTheM
64 std::vector<G4double> G4eDPWAElasticDCS::gTheM
65 std::vector<G4double> G4eDPWAElasticDCS::gTheU
66 std::vector<G4double> G4eDPWAElasticDCS::gTheU
67 // abscissas and weights of an 8 point Gauss-L
68 // for numerical integration on [0,1]
69 const G4double G4eDPWAElasticDCS::gXGL[
70 1.98550718E-02, 1.01666761E-01, 2.37233795E-
71 5.91717321E-01, 7.62766205E-01, 8.98333239E-
72 };
73 const G4double G4eDPWAElasticDCS::gWGL[
74 5.06142681E-02, 1.11190517E-01, 1.56853323E-
75 1.81341892E-01, 1.56853323E-01, 1.11190517E-
76 };
77
78
79 // - iselectron : data for e- (for e+ otherw
80 // - isrestricted : sampling of angular deflec
81 // required (i.e. in case of
82 G4eDPWAElasticDCS::G4eDPWAElasticDCS(G4bool is
83 : fIsRestrictedSamplingRequired(isrestricted),
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
99 if (fSamplingTables[i]) delete fSamplingTa
100 }
101 // clear scp correction data
102 for (std::size_t imc=0; imc<fSCPCPerMatCuts.
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 i
114 void G4eDPWAElasticDCS::InitialiseForZ(std::si
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 fo
124 // should be called only by the master
125 void G4eDPWAElasticDCS::LoadGrid() {
126 G4String fname = FindDirectoryPath() + "grid
127 std::ifstream infile(fname.c_str());
128 if (!infile.is_open()) {
129 G4String msg =
130 " Problem while trying to read " +
131 " G4LEDATA version should be G4EML
132 G4Exception("G4eDPWAElasticDCS::ReadCompre
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)
148 }
149 ++gIndxEnergyLim;
150 // store/set usefull logarithms of the kinet
151 gLogMinEkin = gTheEnergies[0];
152 gInvDelLogEkin = (gNumEnergies-1)/(gTheEnerg
153 // - theta1 in [deg.] (247): we store mu(the
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
160 gTheU1[it] = (theA+1.0)*gTheMus1[it]/(th
161 }
162 // - theta2 in [deg.] (128): we store mu(the
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
168 gTheU2[it] = (theA+1.0)*gTheMus2[it]/(th
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 gNumEnergi
185 const std::size_t hNumEnergries = gNumEner
186 auto v2DHigh = new G4Physics2DVector(gNumT
187 v2DHigh->SetBicubicInterpolation(true);
188 for (std::size_t it=0; it<gNumThetas2; ++i
189 v2DHigh->PutX(it, gTheMus2[it]);
190 }
191 for (std::size_t ie=0; ie<hNumEnergries; +
192 v2DHigh->PutY(ie, gTheEnergies[gIndxEner
193 }
194 std::ostringstream ossh;
195 ossh << FindDirectoryPath() << "dcss/el/dc
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; ++i
200 finh >> dum;
201 for (std::size_t ie=0; ie<hNumEnergries;
202 finh >> dum;
203 v2DHigh->PutValue(it, ie, G4Log(dum*CL
204 }
205 }
206 // load the low energy part:
207 // - with gNumThetas1 theta and gIndxEnerg
208 // for including the firts DCS from the
209 // able to perform interpolation between
210 auto v2DLow = new G4Physics2DVector(gNumTh
211 v2DLow->SetBicubicInterpolation(true);
212 for (std::size_t it=0; it<gNumThetas1; ++i
213 v2DLow->PutX(it, gTheMus1[it]);
214 }
215 for (std::size_t ie=0; ie<gIndxEnergyLim+1
216 v2DLow->PutY(ie, gTheEnergies[ie]);
217 }
218 std::ostringstream ossl;
219 ossl << FindDirectoryPath() << "dcss/el/dc
220 std::istringstream finl(std::ios::in);
221 ReadCompressedFile(ossl.str(), finl);
222 for (std::size_t it=0; it<gNumThetas1; ++i
223 finl >> dum;
224 for (std::size_t ie=0; ie<gIndxEnergyLim
225 finl >> dum;
226 v2DLow->PutValue(it, ie, G4Log(dum*CLH
227 }
228 }
229 // add the +1 part: interpolate the firts
230 std::size_t ix = 0;
231 std::size_t iy = 0;
232 for (std::size_t it=0; it<gNumThetas1; ++i
233 const G4double val = v2DHigh->Value(gThe
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(gNumThetas
242 v2D->SetBicubicInterpolation(true);
243 for (std::size_t it=0; it<gNumThetas2; ++i
244 v2D->PutX(it, gTheMus2[it]);
245 }
246 for (std::size_t ie=0; ie<gNumEnergies; ++
247 v2D->PutY(ie, gTheEnergies[ie]);
248 }
249 std::ostringstream oss;
250 oss << FindDirectoryPath() << "dcss/pos/dc
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; ++i
255 fin >> dum;
256 for (std::size_t ie=0; ie<gNumEnergies;
257 fin >> dum;
258 v2D->PutValue(it, ie, G4Log(dum*CLHEP:
259 }
260 }
261 fDCS[iz]= v2D;
262 }
263 }
264
265
266
267 // Computes the elastic, first and second cros
268 // energy and target atom.
269 // Cross sections are zero ff ekin is below/ab
270 void G4eDPWAElasticDCS::ComputeCSPerAtom(G4int
271 G4dou
272 G4dou
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(thet
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
282 const G4double lekin = std::max(gTheEnergies
283 // if the lower, denser in theta, DCS set sh
284 const G4bool isLowerGrid = (fIsElectron && l
285 const std::vector<G4double>& theMuVector = (
286 const G4Physics2DVector* the2DDCS = (
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)
290 // 0.0 < mumax <= 1.0 for sure here
291 const std::size_t iMuEnd = (mumax == 1.0)
292 // perform numerical integration of the DCS
293 // interval (where mu(theta) = 0.5[1-cos(the
294 std::size_t ix = 0;
295 std::size_t iy = 0;
296 for (std::size_t imu=iMuStart; imu<=iMuEnd;
297 G4double elcsPar = 0.0;
298 G4double tr1csPar = 0.0;
299 G4double tr2csPar = 0.0;
300 const G4double low = (imu==iMuStart) ? mum
301 const G4double del = (imu==iMuEnd) ? mum
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
306 elcsPar += gWGL[igl]*dcs; /
307 tr1csPar += gWGL[igl]*dcs*mu; /
308 tr2csPar += gWGL[igl]*dcs*mu*(1.0-mu); /
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:
321 // NOTE: when Alias is used, sampling on a res
322 // However, Alias makes possible faster
323 // of single scattering model while it's
324 // when restricted interval sampling is
325 // the fIsRestrictedSamplingRequired fla
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; // # da
343 G4double fScreenParA; // the
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 t
353 void G4eDPWAElasticDCS::BuildSmplingTableForZ(
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<OneSamplingTa
359 // read compressed sampling table data
360 std::ostringstream oss;
361 const G4String fname = fIsElectron ? "stable
362 oss << FindDirectoryPath() << fname << "stab
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, !fIsRestrictedSampli
371 // the A screening parameter value used fo
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
400 // interaction (Coulomb scattering) of the pro
401 // fIsElectron) with kinetic energy of exp('le
402 // muber of 'iz'. See the 'SampleCosineThetaRe
403 // a restricted inteval.
404 G4double
405 G4eDPWAElasticDCS::SampleCosineTheta(std::size
406 G4double
407 lekin = std::max(gTheEnergies[0], std::min(g
408 // determine the discrete ekin sampling tabl
409 // - statistical interpolation (i.e. linear
410 const G4double rem = (lekin-gLogMinEkin
411 const auto k = (std::size_t)rem;
412 const std::size_t iekin = (r1 < rem-k) ? k+1
413 // sample the mu(t)=0.5(1-cos(t))
414 const double mu = SampleMu(iz, iekin, r2,
415 return std::max(-1.0, std::min(1.0, 1.0-2.0*
416 }
417
418
419 // samples cos(theta) i.e. cosine of the polar
420 // interaction (Coulomb scattering) of the pro
421 // fIsElectron) with kinetic energy of exp('le
422 // muber of 'iz'.
423 // The cosine theta will be in the [costMin, c
424 // corresponds to a maximum allowed polar scat
425 // costMin corresponds to minimum allowed pola
426 // See the 'SampleCosineTheta' for obtain samp
427 G4double
428 G4eDPWAElasticDCS::SampleCosineThetaRestricted
429
430
431 // costMin corresponds to mu-max while costM
432 lekin = std::max(gTheEnergies[0], std::min(g
433 // determine the discrete ekin sampling tabl
434 // - statistical interpolation (i.e. linear
435 const G4double rem = (lekin-gLogMinEkin
436 const auto k = (size_t)rem;
437 const std::size_t iekin = (r1 < rem-k) ? k :
438 // sample the mu(t)=0.5(1-cos(t))
439 const G4double mu = SampleMu(iz, iekin, r2,
440 return std::max(-1.0, std::min(1.0, 1.0-2.0*
441 }
442
443
444
445 G4double
446 G4eDPWAElasticDCS::SampleMu(std::size_t izet,
447 OneSamplingTable& rtn = (*fSamplingTables[iz
448 // get the lower index of the bin by using t
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[ind
453 // sample value within the selected bin by u
454 const G4double delta = rtn.fCum[indxl + 1] -
455 const G4double aval = r2 * delta;
456
457 const G4double dum1 = (1.0 + rtn.fA[indxl]
458 const G4double dum2 = delta * delta + rtn.f
459 const std::vector<G4double>& theUVect = (fIs
460 const G4double u = theUVect[indxl] + dum
461 // transform back u to mu
462 return rtn.fScreenParA*u/(rtn.fScreenParA+1.
463 }
464
465
466
467 G4double
468 G4eDPWAElasticDCS::FindCumValue(G4double u, co
469 const std::vec
470 const std::size_t iLow = std::distance( uvec
471 const G4double tau = (u-uvect[iLow])/(uv
472 const G4double dum0 = (1.0+stable.fA[iLow
473 const G4double dum1 = 2.0*stable.fB[iLow]
474 const G4double dum2 = 1.0 - std::sqrt(std
475 return std::min(stable.fCum[iLow+1], std::ma
476 }
477
478 // muMin and muMax : no checks on these
479 G4double G4eDPWAElasticDCS::SampleMu(std::size
480 G4double
481 const OneSamplingTable& rtn = (*fSamplingTab
482 const G4double theA = rtn.fScreenParA;
483 //
484 const std::vector<G4double>& theUVect = (fIs
485 const G4double xiMin = (muMin > 0.0) ? Fin
486 const G4double xiMax = (muMax < 1.0) ? Fin
487 //
488 const G4double xi = xiMin+r1*(xiMax-xiM
489 const std::size_t iLow = std::distance( rtn.
490 const G4double delta = rtn.fCum[iLow + 1]
491 const G4double aval = xi - rtn.fCum[iLow]
492
493 const G4double dum1 = (1.0 + rtn.fA[iLow]
494 const G4double dum2 = delta * delta + rtn
495 const G4double u = theUVect[iLow] + du
496 return theA*u/(theA+1.0-u);
497 }
498
499
500
501 // set the DCS data directory path
502 const G4String& G4eDPWAElasticDCS::FindDirecto
503 // check environment variable
504 if (gDataDirectory.empty()) {
505 std::ostringstream ost;
506 ost << G4EmParameters::Instance()->GetDirL
507 gDataDirectory = ost.str();
508 }
509 return gDataDirectory;
510 }
511
512
513
514 // uncompress one data file into the input str
515 void
516 G4eDPWAElasticDCS::ReadCompressedFile(G4String
517 G4String *dataString = nullptr;
518 G4String compfilename(fname+".z");
519 // create input stream with binary mode oper
520 std::ifstream in(compfilename, std::ios::bin
521 if (in.good()) {
522 // get current position in the stream (wa
523 std::size_t fileSize = in.tellg();
524 // set current position being the beginni
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 unco
532 uLongf complen = (uLongf)(fileSize*4);
533 Bytef *uncompdata = new Bytef[complen];
534 while (Z_OK!=uncompress(uncompdata, &comp
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
543 dataString = new G4String((char*)uncompda
544 // delete the uncompressed data buffer
545 delete [] uncompdata;
546 } else {
547 G4String msg =
548 " Problem while trying to read " +
549 " G4LEDATA version should be G4EML
550 G4Exception("G4eDPWAElasticDCS::ReadCompre
551 FatalException,msg.c_str());
552 return;
553 }
554 // create the input string stream from the d
555 if (dataString) {
556 iss.str(*dataString);
557 in.close();
558 delete dataString;
559 }
560 }
561
562
563
564
565 G4double
566 G4eDPWAElasticDCS::ComputeScatteringPowerCorre
567
568 const G4int imc = matcut->GetIndex();
569 G4double corFactor = 1.0;
570 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<
571 return corFactor;
572 }
573 // get the scattering power correction facto
574 const G4double lekin = G4Log(ekin);
575 G4double remaining = (lekin-fSCPCPerMatCut
576 std::size_t lindx = (G4int)remaining;
577 remaining -= lindx;
578 std::size_t imax = fSCPCPerMatCuts[imc]-
579 if (lindx>=imax) {
580 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i
581 } else {
582 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l
583 }
584 return corFactor;
585 }
586
587
588 void G4eDPWAElasticDCS::InitSCPCorrection(G4do
589 G4do
590 // get the material-cuts table
591 G4ProductionCutsTable *thePCTable = G4Produc
592 std::size_t numMatCuts = thePCTab
593 // clear container if any
594 for (std::size_t imc=0; imc<fSCPCPerMatCuts.
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
603 fSCPCPerMatCuts.resize(numMatCuts,nullptr);
604 // loop over the material-cuts and create sc
605 for (G4int imc=0; imc<(G4int)numMatCuts; ++i
606 const G4MaterialCutsCouple *matCut = theP
607 const G4Material* mat = matCut->GetMateri
608 // get e- production cut in the current ma
609 const G4double ecut = (*(thePCTable->Ge
610 const G4double limit = fIsElectron ? 2.0
611 const G4double min = std::max(limit,lo
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
620 numEbins = std::max(numEbins
621 const G4double lmin = G4Log(min);
622 const G4double ldel = G4Log(max/min)/(n
623 fSCPCPerMatCuts[imc] = new SCPCorrection
624 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi
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 param
630 G4double moliereBc = 0.0;
631 G4double moliereXc2 = 0.0;
632 ComputeMParams(mat, moliereBc, moliereXc2)
633 // compute scattering power correction ove
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
638 if (ie>0) {
639 const G4double tau = ekin/CLHEP::
640 const G4double tauCut = ecut/CLHEP::
641 // Moliere's screening parameter
642 const G4double A = moliereXc2/(
643 const G4double gr = (1.+2.*A)*G4
644 const G4double dum0 = (tau+2.)/(ta
645 const G4double dum1 = tau+1.;
646 G4double gm = G4Log(0.5*tau/tauCut)
647 - 0.25*(tau+2.)*( tau+
648 G4Log((tau+4.)*(tau-
649 + 0.5*(tau-2*tauCut)*(
650 if (gm<gr) {
651 gm = gm/gr;
652 } else {
653 gm = 1.;
654 }
655 const G4double z0 = matCut->GetMateri
656 scpCorr = 1.-gm*z0/(z0*(z0+1.));
657 }
658 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo
659 }
660 }
661 }
662
663
664 // compute material dependent Moliere MSC para
665 void G4eDPWAElasticDCS::ComputeMParams(const G
666 G4doubl
667 const G4double const1 = 7821.6; // [
668 const G4double const2 = 0.1569; // [
669 const G4double finstrc2 = 5.325135453E-5; /
670 // G4double xi = 1.0;
671 const G4ElementVector* theElemVect = ma
672 const std::size_t numelems = ma
673 //
674 const G4double* theNbAtomsPerVolVect = ma
675 G4double theTotNbAtomsPerVol = ma
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
683 const G4double zet = (*theElemVect)[iele
684 const G4double iwa = (*theElemVect)[iele
685 const G4double ipz = theNbAtomsPerVolVec
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*finstr
690 sa += ipz*iwa;
691 }
692 const G4double density = mat->GetDensity()*
693 //
694 G4double z0 = (0.0 == sa) ? 0.0 : zs/sa;
695 G4double z1 = (0.0 == zs) ? 0.0 : (ze - zx)
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
700 theBc *= 1.0/CLHEP::cm;
701 theXc2 *= CLHEP::MeV*CLHEP::MeV/CLHEP::cm;
702 }
703
704