| 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 implementation file
30 //
31 // File name: G4GoudsmitSaundersonTable
32 //
33 // Author: Mihaly Novak / (Omrane Kadri
34 //
35 // Creation date: 20.02.2009
36 //
37 // Class description:
38 // Class to handle multiple scattering angul
39 // using Kawrakow-Bielajew Goudsmit-Saunders
40 // Rutherford DCS for elastic scattering of
41 // class is used by G4GoudsmitSaundersonMscM
42 // deflection of electrons/positrons after t
43 //
44 // Modifications:
45 // 04.03.2009 V.Ivanchenko cleanup and format
46 // 26.08.2009 O.Kadri: avoiding unuseful calcu
47 // finding parameter error
48 // 08.02.2010 O.Kadri: reduce delared variable
49 // in secant method
50 // 26.03.2010 O.Kadri: minimum of used arrays
51 // finding method the erro
52 // 12.05.2010 O.Kadri: changing of sqrt((b-a)*
53 // 18.05.2015 M. Novak This class has been com
54 // class name was kept; class descr
55 // A new version of Kawrakow-Bielaj
56 // based on the screened Rutherford
57 // electrons/positrons has been int
58 // angular distributions over a 2D
59 // and the CDFs are now stored in a
60 // together with the corresponding
61 // The new version is several times
62 // compared to the earlier version
63 // that use these data has been als
64 // 28.04.2017 M. Novak: New representation of
65 // significantly reduced data size.
66 // 23.08.2017 M. Novak: Added funtionality to
67 // base GS angular distributions an
68 // parameter, first and second mome
69 // activated in the GS-MSC model.
70 //
71 // References:
72 // [1] A.F.Bielajew, NIMB, 111 (1996) 195-20
73 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19
74 //
75 // -------------------------------------------
76
77 #include "G4GoudsmitSaundersonTable.hh"
78
79
80 #include "G4PhysicalConstants.hh"
81 #include "Randomize.hh"
82 #include "G4Log.hh"
83 #include "G4Exp.hh"
84
85 #include "G4GSMottCorrection.hh"
86 #include "G4MaterialTable.hh"
87 #include "G4Material.hh"
88 #include "G4MaterialCutsCouple.hh"
89 #include "G4ProductionCutsTable.hh"
90 #include "G4EmParameters.hh"
91
92 #include "G4String.hh"
93
94 #include <fstream>
95 #include <cstdlib>
96 #include <cmath>
97
98 #include <iostream>
99 #include <iomanip>
100
101 // perecomputed GS angular distributions, base
102 // are the same for e- and e+ so make sure we
103 G4bool G4GoudsmitSaundersonTable::gIsInitialis
104 //
105 std::vector<G4GoudsmitSaundersonTable::GSMSCAn
106 std::vector<G4GoudsmitSaundersonTable::GSMSCAn
107 //
108 std::vector<double> G4GoudsmitSaundersonTable:
109 std::vector<double> G4GoudsmitSaundersonTable:
110
111
112 G4GoudsmitSaundersonTable::G4GoudsmitSaunderso
113 fIsElectron = iselectron;
114 // set initial values: final values will be
115 fLogLambda0 = 0.; // wi
116 fLogDeltaLambda = 0.; // wi
117 fInvLogDeltaLambda = 0.; // wi
118 fInvDeltaQ1 = 0.; // wi
119 fDeltaQ2 = 0.; // wi
120 fInvDeltaQ2 = 0.; // wi
121 //
122 fLowEnergyLimit = 0.1*CLHEP::keV; // w
123 fHighEnergyLimit = 100.0*CLHEP::MeV; // w
124 //
125 fIsMottCorrection = false; // w
126 fIsPWACorrection = false; // w
127 fMottCorrection = nullptr;
128 //
129 fNumSPCEbinPerDec = 3;
130 }
131
132 G4GoudsmitSaundersonTable::~G4GoudsmitSaunders
133 for (std::size_t i=0; i<gGSMSCAngularDistrib
134 if (gGSMSCAngularDistributions1[i]) {
135 delete [] gGSMSCAngularDistributions1[i]
136 delete [] gGSMSCAngularDistributions1[i]
137 delete [] gGSMSCAngularDistributions1[i]
138 delete gGSMSCAngularDistributions1[i];
139 }
140 }
141 gGSMSCAngularDistributions1.clear();
142 for (std::size_t i=0; i<gGSMSCAngularDistrib
143 if (gGSMSCAngularDistributions2[i]) {
144 delete [] gGSMSCAngularDistributions2[i]
145 delete [] gGSMSCAngularDistributions2[i]
146 delete [] gGSMSCAngularDistributions2[i]
147 delete gGSMSCAngularDistributions2[i];
148 }
149 }
150 gGSMSCAngularDistributions2.clear();
151 if (fMottCorrection) {
152 delete fMottCorrection;
153 fMottCorrection = nullptr;
154 }
155 // clear scp correction data
156 for (std::size_t imc=0; imc<fSCPCPerMatCuts.
157 if (fSCPCPerMatCuts[imc]) {
158 fSCPCPerMatCuts[imc]->fVSCPC.clear();
159 delete fSCPCPerMatCuts[imc];
160 }
161 }
162 fSCPCPerMatCuts.clear();
163 //
164 gIsInitialised = false;
165 }
166
167 void G4GoudsmitSaundersonTable::Initialise(G4d
168 fLowEnergyLimit = lownergylimit;
169 fHighEnergyLimit = highenergylimit;
170 G4double lLambdaMin = G4Log(gLAMBMIN);
171 G4double lLambdaMax = G4Log(gLAMBMAX);
172 fLogLambda0 = lLambdaMin;
173 fLogDeltaLambda = (lLambdaMax-lLambdaMin
174 fInvLogDeltaLambda = 1./fLogDeltaLambda;
175 fInvDeltaQ1 = 1./((gQMAX1-gQMIN1)/(g
176 fDeltaQ2 = (gQMAX2-gQMIN2)/(gQNUM
177 fInvDeltaQ2 = 1./fDeltaQ2;
178 // load precomputed angular distributions an
179 // these are particle independet => they go
180 if (!gIsInitialised) {
181 // load pre-computed GS angular distributi
182 LoadMSCData();
183 gIsInitialised = true;
184 }
185 InitMoliereMSCParams();
186 // Mott-correction: particle(e- or e+) depen
187 if (fIsMottCorrection) {
188 if (!fMottCorrection) {
189 fMottCorrection = new G4GSMottCorrection
190 }
191 fMottCorrection->Initialise();
192 }
193 // init scattering power correction data; us
194 // (Moliere's parameters must be initialised
195 if (fMottCorrection) {
196 InitSCPCorrection();
197 }
198 }
199
200
201 // samplig multiple scattering angles cos(thet
202 // - including no-scattering, single, "few" s
203 // - Mott-correction will be included if it w
204 // lambdaval : s/lambda_el
205 // qval : s/lambda_el G1
206 // scra : screening parameter
207 // cost : will be the smapled cos(theta)
208 // sint : will be the smapled sin(theta)
209 // lekin : logarithm of the current kineti
210 // beta2 : the corresponding beta square
211 // matindx : index of the current material
212 // returns true if it was msc
213 G4bool G4GoudsmitSaundersonTable::Sampling(G4d
214 G4d
215 GSM
216 G4d
217 G4double rand0 = G4UniformRand();
218 G4double expn = G4Exp(-lambdaval);
219 //
220 // no scattering case
221 if (rand0<expn) {
222 cost = 1.0;
223 sint = 0.0;
224 return false;
225 }
226 //
227 // single scattering case : sample from the
228 // - Mott-correction will be included if it
229 if (rand0<(1.+lambdaval)*expn) {
230 // cost is sampled in SingleScattering()
231 cost = SingleScattering(lambdaval, scra, l
232 // add protections
233 if (cost<-1.0) cost = -1.0;
234 if (cost>1.0) cost = 1.0;
235 // compute sin(theta) from the sampled cos
236 G4double dum0 = 1.-cost;
237 sint = std::sqrt(dum0*(2.0-dum0));
238 return false;
239 }
240 //
241 // handle this case:
242 // -lambdaval < 1 i.e. mean #elastic ev
243 // the currently sampled case is not 0
244 // lambdaval (that we have precomputed
245 // stored in a form of equally probabe
246 // interp. parameters) is 1.]
247 // -probability of having n elastic eve
248 // lambdaval parameter.
249 // -the max. probability (when lambdava
250 // elastic events is 0.2642411 and the
251 // events decays rapidly with n. So se
252 // -sampling of this cases is done in a
253 // where the current #elastic event is
254 if (lambdaval<1.0) {
255 G4double prob, cumprob;
256 prob = cumprob = expn;
257 G4double curcost,cursint;
258 // init cos(theta) and sin(theta) to the z
259 cost = 1.0;
260 sint = 0.0;
261 for (G4int iel=1; iel<10; ++iel) {
262 // prob of having iel scattering from Po
263 prob *= lambdaval/(G4double)iel;
264 cumprob += prob;
265 //
266 //sample cos(theta) from the singe scatt
267 // - Mott-correction will be included if
268 curcost = SingleScattering(lambdav
269 G4double dum0 = 1.-curcost;
270 cursint = dum0*(2.0-dum0); // sin^
271 //
272 // if we got current deflection that is
273 // then update cos(theta) sin(theta)
274 if (cursint>1.0e-20) {
275 cursint = std::sqrt(cursint);
276 G4double curphi = CLHEP::twopi*G4Unifo
277 cost = cost*curcost-sint*cu
278 sint = std::sqrt(std::max(0
279 }
280 //
281 // check if we have done enough scatteri
282 if (rand0<cumprob) {
283 return false;
284 }
285 }
286 // if reached the max iter i.e. 10
287 return false;
288 }
289 //
290 // multiple scattering case with lambdavalue
291 // - use the precomputed and transformed G
292 // distributions to sample cos(theta)
293 // - Mott-correction will be included if i
294 cost = SampleCosTheta(lambdaval, qval, scra,
295 // add protections
296 if (cost<-1.0) cost = -1.0;
297 if (cost> 1.0) cost = 1.0;
298 // compute cos(theta) and sin(theta) from th
299 G4double dum0 = 1.0-cost;
300 sint = std::sqrt(dum0*(2.0-dum0));
301 // return true if it was msc
302 return true;
303 }
304
305
306 G4double G4GoudsmitSaundersonTable::SampleCosT
307
308
309
310 G4double cost = 1.;
311 // determine the base GS angular distributio
312 if (isfirst) {
313 *gsDtr = GetGSAngularDtr(scra, lambdaval,
314 }
315 // sample cost from the GS angular distribut
316 cost = SampleGSSRCosTheta(*gsDtr, transfPar)
317 // Mott-correction if it was requested by th
318 if (fIsMottCorrection && *gsDtr) { // no
319 static const G4int nlooplim = 1000;
320 G4int nloop = 0 ; // rejection loop
321 // G4int ekindx = -1; // evaluate only
322 // G4int deltindx = -1 ; // evaluate onl
323 G4double val = fMottCorrection->GetMo
324 while (G4UniformRand()>val && ++nloop<nloo
325 // sampling cos(theta)
326 cost = SampleGSSRCosTheta(*gsDtr, transf
327 val = fMottCorrection->GetMottRejection
328 };
329 }
330 return cost;
331 }
332
333
334 // returns with cost sampled from the GS angul
335 G4double G4GoudsmitSaundersonTable::SampleGSSR
336 // check if isotropic theta (i.e. cost is un
337 if (!gsDtr) {
338 return 1.-2.0*G4UniformRand();
339 }
340 //
341 // sampling form the selected distribution
342 G4double ndatm1 = gsDtr->fNumData-1.;
343 G4double delta = 1.0/ndatm1;
344 // determine lower cumulative bin inidex
345 G4double rndm = G4UniformRand();
346 G4int indxl = rndm*ndatm1;
347 G4double aval = rndm-indxl*delta;
348 G4double dum0 = delta*aval;
349
350 G4double dum1 = (1.0+gsDtr->fParamA[indxl]
351 G4double dum2 = delta*delta + gsDtr->fPara
352 G4double sample = gsDtr->fUValues[indxl] +
353 // transform back u to cos(theta) :
354 // this is the sampled cos(theta) = (2.0*par
355 return 1.-(2.0*transfpar*sample)/(1.0-sample
356 }
357
358
359 // determine the GS angular distribution we ne
360 G4GoudsmitSaundersonTable::GSMSCAngularDtr* G4
361
362 GSMSCAngularDtr *dtr = nullptr;
363 G4bool first = false;
364 // isotropic cost above gQMAX2 (i.e. dtr sta
365 if (qval<gQMAX2) {
366 G4int lamIndx = -1; // lambda value in
367 G4int qIndx = -1; // lambda value in
368 // init to second grid Q values
369 G4int numQVal = gQNUM2;
370 G4double minQVal = gQMIN2;
371 G4double invDelQ = fInvDeltaQ2;
372 G4double pIndxH = 0.; // probability of
373 // check if first or second grid needs to
374 if (qval<gQMIN2) { // first grid
375 first = true;
376 // protect against qval<gQMIN1
377 if (qval<gQMIN1) {
378 qval = gQMIN1;
379 qIndx = 0;
380 //pIndxH = 0.;
381 }
382 // set to first grid Q values
383 numQVal = gQNUM1;
384 minQVal = gQMIN1;
385 invDelQ = fInvDeltaQ1;
386 }
387 // make sure that lambda = s/lambda_el is
388 // lambda<gLAMBMIN=1 is already handeled b
389 if (lambdaval>=gLAMBMAX) {
390 lambdaval = gLAMBMAX-1.e-8;
391 lamIndx = gLAMBNUM-1;
392 }
393 G4double lLambda = G4Log(lambdaval);
394 //
395 // determine lower lambda (=s/lambda_el) i
396 if (lamIndx<0) {
397 pIndxH = (lLambda-fLogLambda0)*fInvLogD
398 lamIndx = (G4int)(pIndxH); // low
399 pIndxH = pIndxH-lamIndx; // proba
400 if (G4UniformRand()<pIndxH) {
401 ++lamIndx;
402 }
403 }
404 //
405 // determine lower Q (=s/lambda_el G1) ind
406 if (qIndx<0) {
407 pIndxH = (qval-minQVal)*invDelQ;
408 qIndx = (G4int)(pIndxH); // low
409 pIndxH = pIndxH-qIndx;
410 if (G4UniformRand()<pIndxH) {
411 ++qIndx;
412 }
413 }
414 // set indx
415 G4int indx = lamIndx*numQVal+qIndx;
416 if (first) {
417 dtr = gGSMSCAngularDistributions1[indx];
418 } else {
419 dtr = gGSMSCAngularDistributions2[indx];
420 }
421 // dtr might be nullptr that indicates iso
422 // - if the selected lamIndx, qIndx corres
423 // G1 should always be < 1 and if G1 is
424 //
425 // compute the transformation parameter
426 if (lambdaval>10.0) {
427 transfpar = 0.5*(-2.77164+lLambda*( 2.94
428 } else {
429 transfpar = 0.5*(1.347+lLambda*(0.209364
430 }
431 transfpar *= (lambdaval+4.0)*scra;
432 }
433 // return with the selected GS angular distr
434 return dtr;
435 }
436
437
438 void G4GoudsmitSaundersonTable::LoadMSCData()
439 gGSMSCAngularDistributions1.resize(gLAMBNUM*
440 const G4String str1 = G4EmParameters::Instan
441 for (G4int il=0; il<gLAMBNUM; ++il) {
442 G4String fname = str1 + std::to_string(il)
443 std::ifstream infile(fname,std::ios::in);
444 if (!infile.is_open()) {
445 G4String msgc = "Cannot open file: " + f
446 G4Exception("G4GoudsmitSaundersonTable::
447 FatalException, msgc.c_str());
448 return;
449 }
450 for (G4int iq=0; iq<gQNUM1; ++iq) {
451 auto gsd = new GSMSCAngularDtr();
452 infile >> gsd->fNumData;
453 gsd->fUValues = new G4double[gsd->fNumDa
454 gsd->fParamA = new G4double[gsd->fNumDa
455 gsd->fParamB = new G4double[gsd->fNumDa
456 G4double ddummy;
457 infile >> ddummy; infile >> ddummy;
458 for (G4int i=0; i<gsd->fNumData; ++i) {
459 infile >> gsd->fUValues[i];
460 infile >> gsd->fParamA[i];
461 infile >> gsd->fParamB[i];
462 }
463 gGSMSCAngularDistributions1[il*gQNUM1+iq
464 }
465 infile.close();
466 }
467 //
468 // second grid
469 gGSMSCAngularDistributions2.resize(gLAMBNUM*
470 const G4String str2 = G4EmParameters::Instan
471 for (G4int il=0; il<gLAMBNUM; ++il) {
472 G4String fname = str2 + std::to_string(il)
473 std::ifstream infile(fname,std::ios::in);
474 if (!infile.is_open()) {
475 G4String msgc = "Cannot open file: " + f
476 G4Exception("G4GoudsmitSaundersonTable::
477 FatalException, msgc.c_str());
478 return;
479 }
480 for (G4int iq=0; iq<gQNUM2; ++iq) {
481 G4int numData;
482 infile >> numData;
483 if (numData>1) {
484 auto gsd = new GSMSCAngularDtr();
485 gsd->fNumData = numData;
486 gsd->fUValues = new G4double[gsd->fNum
487 gsd->fParamA = new G4double[gsd->fNum
488 gsd->fParamB = new G4double[gsd->fNum
489 double ddummy;
490 infile >> ddummy; infile >> ddummy;
491 for (G4int i=0; i<gsd->fNumData; ++i)
492 infile >> gsd->fUValues[i];
493 infile >> gsd->fParamA[i];
494 infile >> gsd->fParamB[i];
495 }
496 gGSMSCAngularDistributions2[il*gQNUM2+
497 } else {
498 gGSMSCAngularDistributions2[il*gQNUM2+
499 }
500 }
501 infile.close();
502 }
503 }
504
505 // samples cost in single scattering based on
506 // (with Mott-correction if it was requested)
507 G4double G4GoudsmitSaundersonTable::SingleScat
508
509
510 G4double rand1 = G4UniformRand();
511 // sample cost from the Screened-Rutherford
512 G4double cost = 1.-2.0*scra*rand1/(1.0-rand
513 // Mott-correction if it was requested by th
514 if (fIsMottCorrection) {
515 static const G4int nlooplim = 1000; // re
516 G4int nloop = 0 ; // loop counter
517 G4int ekindx = -1 ; // evaluate only
518 G4int deltindx = 0 ; // single scatter
519 G4double q1 = 0.; // not used when
520 // computing Mott rejection function value
521 G4double val = fMottCorrection->GetMo
522
523 while (G4UniformRand()>val && ++nloop<nloo
524 // sampling cos(theta) from the Screened
525 rand1 = G4UniformRand();
526 cost = 1.-2.0*scra*rand1/(1.0-rand1+scr
527 // computing Mott rejection function val
528 val = fMottCorrection->GetMottRejectio
529
530 };
531 }
532 return cost;
533 }
534
535
536 void G4GoudsmitSaundersonTable::GetMottCorrec
537
538
539 if (fIsMottCorrection) {
540 fMottCorrection->GetMottCorrectionFactors(
541 }
542 }
543
544
545 // compute material dependent Moliere MSC para
546 void G4GoudsmitSaundersonTable::InitMoliereMSC
547 const G4double const1 = 7821.6; // [
548 const G4double const2 = 0.1569; // [
549 const G4double finstrc2 = 5.325135453E-5; /
550
551 G4MaterialTable* theMaterialTable = G4Mater
552 // get number of materials in the table
553 std::size_t numMaterials = theMaterialTable
554 // make sure that we have long enough vecto
555 if(gMoliereBc.size()<numMaterials) {
556 gMoliereBc.resize(numMaterials);
557 gMoliereXc2.resize(numMaterials);
558 }
559 G4double xi = 1.0;
560 G4int maxZ = 200;
561 if (fIsMottCorrection || fIsPWACorrection)
562 // xi = 1.0; <= always set to 1 from n
563 maxZ = G4GSMottCorrection::GetMaxZet();
564 }
565 //
566 for (std::size_t imat=0; imat<numMaterials;
567 const G4Material* theMaterial =
568 const G4ElementVector* theElemVect =
569 const G4int numelems =
570 //
571 const G4double* theNbAtomsPerVolVe
572 G4double theTotNbAtomsPerVo
573 //
574 G4double zs = 0.0;
575 G4double zx = 0.0;
576 G4double ze = 0.0;
577 G4double sa = 0.0;
578 //
579 for(G4int ielem = 0; ielem < numelems; ie
580 G4double zet = (*theElemVect)[ielem]->G
581 if (zet>maxZ) {
582 zet = (G4double)maxZ;
583 }
584 G4double iwa = (*theElemVect)[ielem]->
585 G4double ipz = theNbAtomsPerVolVect[ie
586 G4double dum = ipz*zet*(zet+xi);
587 zs += dum;
588 ze += dum*(-2.0/3.0)*G4Log(ze
589 zx += dum*G4Log(1.0+3.34*fins
590 sa += ipz*iwa;
591 }
592 G4double density = theMaterial->GetDensit
593 //
594 gMoliereBc[theMaterial->GetIndex()] = co
595 gMoliereXc2[theMaterial->GetIndex()] = co
596 // change to Geant4 internal units of 1/l
597 gMoliereBc[theMaterial->GetIndex()] *= 1
598 gMoliereXc2[theMaterial->GetIndex()] *= C
599 }
600 }
601
602
603 // this method is temporary, will be removed/r
604 G4double G4GoudsmitSaundersonTable::ComputeSca
605 G4int imc = matcut->GetIndex();
606 G4double corFactor = 1.0;
607 if (!(fSCPCPerMatCuts[imc]->fIsUse) || ekin<
608 return corFactor;
609 }
610 // get the scattering power correction facto
611 G4double lekin = G4Log(ekin);
612 G4double remaining = (lekin-fSCPCPerMatCuts
613 std::size_t lindx = (std::size_t)remaining
614 remaining -= lindx;
615 std::size_t imax = fSCPCPerMatCuts[imc]->
616 if (lindx>=imax) {
617 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[i
618 } else {
619 corFactor = fSCPCPerMatCuts[imc]->fVSCPC[l
620 }
621 return corFactor;
622 }
623
624
625 void G4GoudsmitSaundersonTable::InitSCPCorrect
626 // get the material-cuts table
627 G4ProductionCutsTable *thePCTable = G4Produc
628 std::size_t numMatCuts = thePCTab
629 // clear container if any
630 for (std::size_t imc=0; imc<fSCPCPerMatCuts.
631 if (fSCPCPerMatCuts[imc]) {
632 fSCPCPerMatCuts[imc]->fVSCPC.clear();
633 delete fSCPCPerMatCuts[imc];
634 fSCPCPerMatCuts[imc] = nullptr;
635 }
636 }
637 //
638 // set size of the container and create the
639 fSCPCPerMatCuts.resize(numMatCuts,nullptr);
640 // loop over the material-cuts and create sc
641 for (G4int imc=0; imc<(G4int)numMatCuts; ++i
642 const G4MaterialCutsCouple *matCut = theP
643 // get e- production cut in the current ma
644 G4double limit;
645 G4double ecut;
646 if (fIsElectron) {
647 ecut = (*(thePCTable->GetEnergyCutsVect
648 limit = 2.*ecut;
649 } else {
650 ecut = (*(thePCTable->GetEnergyCutsVect
651 limit = ecut;
652 }
653 G4double min = std::max(limit,fLowEnergyLi
654 G4double max = fHighEnergyLimit;
655 if (min>=max) {
656 fSCPCPerMatCuts[imc] = new SCPCorrection
657 fSCPCPerMatCuts[imc]->fIsUse = false;
658 fSCPCPerMatCuts[imc]->fPrCut = min;
659 continue;
660 }
661 G4int numEbins = fNumSPCEbinPerDec*G
662 numEbins = std::max(numEbins,3
663 G4double lmin = G4Log(min);
664 G4double ldel = G4Log(max/min)/(num
665 fSCPCPerMatCuts[imc] = new SCPCorrection()
666 fSCPCPerMatCuts[imc]->fVSCPC.resize(numEbi
667 fSCPCPerMatCuts[imc]->fIsUse = true;
668 fSCPCPerMatCuts[imc]->fPrCut = min;
669 fSCPCPerMatCuts[imc]->fLEmin = lmin;
670 fSCPCPerMatCuts[imc]->fILDel = 1./ldel;
671 for (G4int ie=0; ie<numEbins; ++ie) {
672 G4double ekin = G4Exp(lmin+ie*ldel);
673 G4double scpCorr = 1.0;
674 // compute correction factor: I.Kawrakow
675 if (ie>0) {
676 G4double tau = ekin/CLHEP::electr
677 G4double tauCut = ecut/CLHEP::electr
678 // Moliere's screening parameter
679 G4int matindx = (G4int)matCut->Get
680 G4double A = GetMoliereXc2(mati
681 G4double gr = (1.+2.*A)*G4Log(1.
682 G4double dum0 = (tau+2.)/(tau+1.);
683 G4double dum1 = tau+1.;
684 G4double gm = G4Log(0.5*tau/tauC
685 - 0.25*(tau+2.)*(
686 G4Log((tau+4.)*(
687 + 0.5*(tau-2*tauCu
688 if (gm<gr) {
689 gm = gm/gr;
690 } else {
691 gm = 1.;
692 }
693 G4double z0 = matCut->GetMaterial()->
694 scpCorr = 1.-gm*z0/(z0*(z0+1.));
695 }
696 fSCPCPerMatCuts[imc]->fVSCPC[ie] = scpCo
697 }
698 }
699 }
700