| 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
30 //
31 // File name: G4EmCorrections
32 //
33 // Author: Vladimir Ivanchenko
34 //
35 // Creation date: 13.01.2005
36 //
37 // Modifications:
38 // 05.05.2005 V.Ivanchenko Fix misprint in Mot
39 // 26.11.2005 V.Ivanchenko Fix effective charg
40 // 28.04.2006 V.Ivanchenko General cleanup, ad
41 // 13.05.2006 V.Ivanchenko Add corrections for
42 // 08.05.2007 V.Ivanchenko Use G4IonTable for
43 // division by zero
44 // 29.02.2008 V.Ivanchenko use expantions for
45 // 21.04.2008 Updated computations for ions (V
46 // 20.05.2008 Removed Finite Size correction (
47 // 19.04.2012 Fix reproducibility problem (A.R
48 //
49 //
50 // Class Description:
51 //
52 // This class provides calculation of EM corre
53 //
54
55 // -------------------------------------------
56 //
57
58 #include "G4EmCorrections.hh"
59 #include "G4PhysicalConstants.hh"
60 #include "G4SystemOfUnits.hh"
61 #include "G4ParticleTable.hh"
62 #include "G4IonTable.hh"
63 #include "G4VEmModel.hh"
64 #include "G4Proton.hh"
65 #include "G4GenericIon.hh"
66 #include "G4PhysicsLogVector.hh"
67 #include "G4ProductionCutsTable.hh"
68 #include "G4MaterialCutsCouple.hh"
69 #include "G4AtomicShells.hh"
70 #include "G4Log.hh"
71 #include "G4Exp.hh"
72 #include "G4Pow.hh"
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 namespace
77 {
78 constexpr G4double inveplus = 1.0/CLHEP::epl
79 constexpr G4double alpha2 = CLHEP::fine_stru
80 }
81 const G4double G4EmCorrections::ZD[11] =
82 {0., 0., 0., 1.72, 2.09, 2.48, 2.82, 3.16,
83 const G4double G4EmCorrections::UK[20] = {1.99
84 2.0817, 2.0945, 2.0
85 2.1197, 2.1246, 2.1
86 2.1310, 2.1310, 2.1
87 const G4double G4EmCorrections::VK[20] = {8.34
88 8.3219, 8.3201, 8.3
89 8.3191, 8.3199, 8.3
90 8.3244, 8.3264, 8.3
91 G4double G4EmCorrections::ZK[] = {0.0};
92 const G4double G4EmCorrections::Eta[29] = {0.0
93 0.03, 0.04, 0.05
94 0.1, 0.15, 0.2,
95 0.5, 0.6, 0.7,
96 1.2, 1.4, 1.5,
97 G4double G4EmCorrections::CK[20][29];
98 G4double G4EmCorrections::CL[26][28];
99 const G4double G4EmCorrections::UL[] = {0.1215
100 1.4379, 1.5032, 1.5
101 1.8036, 1.8543, 1.8
102 1.9508, 1.9696, 1.9
103 2.0001, 2.0039, 2.0
104 G4double G4EmCorrections::VL[] = {0.0};
105 G4double G4EmCorrections::sWmaxBarkas = 10.0;
106
107 G4PhysicsFreeVector* G4EmCorrections::sBarkasC
108 G4PhysicsFreeVector* G4EmCorrections::sThetaK
109 G4PhysicsFreeVector* G4EmCorrections::sThetaL
110
111 G4EmCorrections::G4EmCorrections(G4int verb)
112 : verbose(verb)
113 {
114 eth = 2.0*CLHEP::MeV;
115 eCorrMin = 25.*CLHEP::keV;
116 eCorrMax = 1.*CLHEP::GeV;
117
118 ionTable = G4ParticleTable::GetParticleTable
119 g4calc = G4Pow::GetInstance();
120
121 // fill vectors
122 if (nullptr == sBarkasCorr) {
123 Initialise();
124 isInitializer = true;
125 }
126 }
127
128 //....oooOO0OOooo........oooOO0OOooo........oo
129
130 G4EmCorrections::~G4EmCorrections()
131 {
132 for (G4int i=0; i<nIons; ++i) { delete stopD
133 if (isInitializer) {
134 delete sBarkasCorr;
135 delete sThetaK;
136 delete sThetaL;
137 sBarkasCorr = sThetaK = sThetaL = nullptr;
138 }
139 }
140
141 void G4EmCorrections::SetupKinematics(const G4
142 const G4Material* mat,
143 const G4double kineticEnergy)
144 {
145 if(kineticEnergy != kinEnergy || p != partic
146 particle = p;
147 kinEnergy = kineticEnergy;
148 mass = p->GetPDGMass();
149 tau = kineticEnergy / mass;
150 gamma = 1.0 + tau;
151 bg2 = tau * (tau+2.0);
152 beta2 = bg2/(gamma*gamma);
153 beta = std::sqrt(beta2);
154 ba2 = beta2/alpha2;
155 const G4double ratio = CLHEP::electron_mas
156 tmax = 2.0*CLHEP::electron_mass_c2*bg2
157 /(1. + 2.0*gamma*ratio + ratio*ratio);
158 charge = p->GetPDGCharge()*inveplus;
159 if(charge > 1.5) { charge = effCharge.Effe
160 q2 = charge*charge;
161 }
162 if(mat != material) {
163 material = mat;
164 theElementVector = material->GetElementVec
165 atomDensity = material->GetAtomicNumDensi
166 numberOfElements = (G4int)material->GetNum
167 }
168 }
169
170 //....oooOO0OOooo........oooOO0OOooo........oo
171
172 G4double G4EmCorrections::HighOrderCorrections
173
174
175 {
176 // . Z^3 Barkas effect in the stopping power
177 // J.C Ashley and R.H.Ritchie
178 // Physical review B Vol.5 No.7 1 April 19
179 // and ICRU49 report
180 // valid for kineticEnergy < 0.5 MeV
181 // Other corrections from S.P.Ahlen Rev. M
182
183 SetupKinematics(p, mat, e);
184 if(tau <= 0.0) { return 0.0; }
185
186 const G4double Barkas = BarkasCorrection(p,
187 const G4double Bloch = BlochCorrection(p, m
188 const G4double Mott = MottCorrection(p, mat,
189
190 G4double sum = 2.0*(Barkas + Bloch) + Mott;
191
192 if(verbose > 1) {
193 G4cout << "EmCorrections: E(MeV)= " << e/M
194 << " Bloch= " << Bloch << " Mott= "
195 << " Sum= " << sum << " q2= " << q2
196 G4cout << " ShellCorrection: " << ShellCor
197 << " Kshell= " << KShellCorrection(
198 << " Lshell= " << LShellCorrection(
199 << " " << mat->GetName() << G4end
200 }
201 sum *= material->GetElectronDensity()*q2*CLH
202 return sum;
203 }
204
205 //....oooOO0OOooo........oooOO0OOooo........oo
206
207 G4double G4EmCorrections::IonBarkasCorrection(
208
209
210 {
211 SetupKinematics(p, mat, e);
212 return 2.0*BarkasCorrection(p, mat, e, true)
213 material->GetElectronDensity() * q2 * CLHE
214 }
215
216 //....oooOO0OOooo........oooOO0OOooo........oo
217
218 G4double G4EmCorrections::ComputeIonCorrection
219
220
221 {
222 // . Z^3 Barkas effect in the stopping power
223 // J.C Ashley and R.H.Ritchie
224 // Physical review B Vol.5 No.7 1 April 19
225 // and ICRU49 report
226 // valid for kineticEnergy < 0.5 MeV
227 // Other corrections from S.P.Ahlen Rev. M
228 SetupKinematics(p, mat, e);
229 if(tau <= 0.0) { return 0.0; }
230
231 const G4double Barkas = BarkasCorrection (p,
232 const G4double Bloch = BlochCorrection (p,
233 const G4double Mott = MottCorrection (p, mat
234
235 G4double sum = 2.0*(Barkas*(charge - 1.0)/ch
236
237 if(verbose > 1) {
238 G4cout << "EmCorrections: E(MeV)= " << e/M
239 << " Bloch= " << Bloch << " Mott= "
240 << " Sum= " << sum << G4endl;
241 }
242 sum *= material->GetElectronDensity() * q2 *
243
244 if(verbose > 1) { G4cout << " Sum= " << sum
245 return sum;
246 }
247
248 //....oooOO0OOooo........oooOO0OOooo........oo
249
250 G4double G4EmCorrections::IonHighOrderCorrecti
251
252
253 {
254 // . Z^3 Barkas effect in the stopping power
255 // J.C Ashley and R.H.Ritchie
256 // Physical review B Vol.5 No.7 1 April 19
257 // and ICRU49 report
258 // valid for kineticEnergy < 0.5 MeV
259 // Other corrections from S.P.Ahlen Rev. M
260
261 G4double sum = 0.0;
262
263 if (nullptr != ionHEModel) {
264 G4int Z = G4lrint(p->GetPDGCharge()*invepl
265 Z = std::max(std::min(Z, 99), 1);
266
267 const G4double ethscaled = eth*p->GetPDGMa
268 const G4int ionPDG = p->GetPDGEncoding();
269 auto iter = thcorr.find(ionPDG);
270 if (iter == thcorr.end()) { // Not found:
271 std::vector<G4double> v;
272 for(std::size_t i=0; i<ncouples; ++i){
273 v.push_back(ethscaled*ComputeIonCorrec
274 }
275 thcorr.insert(std::pair< G4int, std::vec
276 }
277 G4double rest = 0.0;
278 iter = thcorr.find(ionPDG);
279 if (iter != thcorr.end()) { rest = (iter->
280
281 sum = ComputeIonCorrections(p,couple->GetM
282
283 if(verbose > 1) {
284 G4cout << " Sum= " << sum << " dSum= " <
285 }
286 }
287 return sum;
288 }
289
290 //....oooOO0OOooo........oooOO0OOooo........oo
291
292 G4double G4EmCorrections::Bethe(const G4Partic
293 const G4Materi
294 const G4double
295 {
296 SetupKinematics(p, mat, e);
297 const G4double eexc = material->GetIonisati
298 const G4double eexc2 = eexc*eexc;
299 return 0.5*G4Log(2.0*electron_mass_c2*bg2*tm
300 }
301
302 //....oooOO0OOooo........oooOO0OOooo........oo
303
304 G4double G4EmCorrections::SpinCorrection(const
305 const
306 const
307 {
308 SetupKinematics(p, mat, e);
309 const G4double dedx = 0.5*tmax/(kinEnergy +
310 return 0.5*dedx*dedx;
311 }
312
313 //....oooOO0OOooo........oooOO0OOooo........oo
314
315 G4double G4EmCorrections:: KShellCorrection(co
316 co
317 co
318 {
319 SetupKinematics(p, mat, e);
320 G4double term = 0.0;
321 for (G4int i = 0; i<numberOfElements; ++i) {
322
323 G4double Z = (*theElementVector)[i]->GetZ(
324 G4int iz = (*theElementVector)[i]->GetZa
325 G4double f = 1.0;
326 G4double Z2= (Z-0.3)*(Z-0.3);
327 if(1 == iz) {
328 f = 0.5;
329 Z2 = 1.0;
330 }
331 const G4double eta = ba2/Z2;
332 const G4double tet = (11 < iz) ? sThetaK->
333 term += f*atomDensity[i]*KShell(tet,eta)/Z
334 }
335
336 term /= material->GetTotNbOfAtomsPerVolume()
337
338 return term;
339 }
340
341 //....oooOO0OOooo........oooOO0OOooo........oo
342
343 G4double G4EmCorrections:: LShellCorrection(co
344 co
345 co
346 {
347 SetupKinematics(p, mat, e);
348 G4double term = 0.0;
349 for (G4int i = 0; i<numberOfElements; ++i) {
350
351 const G4double Z = (*theElementVector)[i]-
352 const G4int iz = (*theElementVector)[i]->G
353 if(2 < iz) {
354 const G4double Zeff = (iz < 10) ? Z - ZD
355 const G4double Z2= Zeff*Zeff;
356 const G4double eta = ba2/Z2;
357 G4double tet = sThetaL->Value(Z);
358 G4int nmax = std::min(4, G4AtomicShells:
359 for (G4int j=1; j<nmax; ++j) {
360 G4int ne = G4AtomicShells::GetNumberOf
361 if (15 >= iz) {
362 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2*
363 0.25*Z2*(1.0 + Z2*alpha2/16.);
364 }
365 //G4cout << " LShell: j= " << j << " n
366 // << " ThetaL= " << tet << G4en
367 term += 0.125*ne*atomDensity[i]*LShell
368 }
369 }
370 }
371
372 term /= material->GetTotNbOfAtomsPerVolume()
373
374 return term;
375 }
376
377 //....oooOO0OOooo........oooOO0OOooo........oo
378
379 G4double G4EmCorrections::KShell(const G4doubl
380 {
381 G4double corr = 0.0;
382
383 static const G4double TheK[20] =
384 {0.64, 0.65, 0.66, 0.68, 0.70, 0.72, 0.74,
385 0.80, 0.82, 0.84, 0.85, 0.86, 0.88, 0.90,
386
387
388 G4double x = tet;
389 G4int itet = 0;
390 G4int ieta = 0;
391 if(tet < TheK[0]) {
392 x = TheK[0];
393 } else if(tet > TheK[nK-1]) {
394 x = TheK[nK-1];
395 itet = nK-2;
396 } else {
397 itet = Index(x, TheK, nK);
398 }
399 // asymptotic case
400 if(eta >= Eta[nEtaK-1]) {
401 corr =
402 (Value(x, TheK[itet], TheK[itet+1], UK[i
403 Value(x, TheK[itet], TheK[itet+1], VK[i
404 Value(x, TheK[itet], TheK[itet+1], ZK[i
405 } else {
406 G4double y = eta;
407 if(eta < Eta[0]) {
408 y = Eta[0];
409 } else {
410 ieta = Index(y, Eta, nEtaK);
411 }
412 corr = Value2(x, y, TheK[itet], TheK[itet+
413 CK[itet][ieta], CK[itet+1][i
414 CK[itet][ieta+1], CK[itet+1]
415 //G4cout << " x= " <<x<<" y= "<<y<<" tet
416 // <<" "<< TheK[itet+1]<<" eta=
417 // <<" CK= " << CK[itet][ieta]<<
418 // <<" "<< CK[itet][ieta+1]<<" "
419 }
420 //G4cout << "Kshell: tet= " << tet << " eta
421 // << " itet= " << itet << " ieta=
422 return corr;
423 }
424
425 //....oooOO0OOooo........oooOO0OOooo........oo
426
427 G4double G4EmCorrections::LShell(const G4doubl
428 {
429 G4double corr = 0.0;
430
431 static const G4double TheL[26] =
432 {0.24, 0.26, 0.28, 0.30, 0.32, 0.34, 0.3
433 0.42, 0.44, 0.45, 0.46, 0.48, 0.50, 0.5
434 0.58, 0.60, 0.62, 0.64, 0.65, 0.66};
435
436 G4double x = tet;
437 G4int itet = 0;
438 G4int ieta = 0;
439 if(tet < TheL[0]) {
440 x = TheL[0];
441 } else if(tet > TheL[nL-1]) {
442 x = TheL[nL-1];
443 itet = nL-2;
444 } else {
445 itet = Index(x, TheL, nL);
446 }
447
448 // asymptotic case
449 if(eta >= Eta[nEtaL-1]) {
450 corr = (Value(x, TheL[itet], TheL[itet+1],
451 + Value(x, TheL[itet], TheL[itet
452 )/eta;
453 } else {
454 G4double y = eta;
455 if(eta < Eta[0]) {
456 y = Eta[0];
457 } else {
458 ieta = Index(y, Eta, nEtaL);
459 }
460 corr = Value2(x, y, TheL[itet], TheL[itet+
461 CL[itet][ieta], CL[itet+1][i
462 CL[itet][ieta+1], CL[itet+1]
463 //G4cout << " x= " <<x<<" y= "<<y<<" tet
464 // <<" "<< TheL[itet+1]<<" eta=
465 // <<" CL= " << CL[itet][ieta]<<
466 // <<" "<< CL[itet][ieta+1]<<" "
467 }
468 //G4cout<<"Lshell: tet= "<<tet<<" eta= "<<e
469 // <<" ieta= "<<ieta<<" Corr= "<<cor
470 return corr;
471 }
472
473 //....oooOO0OOooo........oooOO0OOooo........oo
474
475 G4double G4EmCorrections::ShellCorrectionSTD(c
476 c
477 c
478 {
479 SetupKinematics(p, mat, e);
480 G4double taulim= 8.0*MeV/mass;
481 G4double bg2lim= taulim * (taulim+2.0);
482
483 G4double* shellCorrectionVector =
484 material->GetIonisation()->GetShel
485 G4double sh = 0.0;
486 G4double x = 1.0;
487 G4double taul = material->GetIonisation()->
488
489 if ( bg2 >= bg2lim ) {
490 for (G4int k=0; k<3; ++k) {
491 x *= bg2 ;
492 sh += shellCorrectionVector[k]/x;
493 }
494
495 } else {
496 for (G4int k=0; k<3; ++k) {
497 x *= bg2lim ;
498 sh += shellCorrectionVector[k]/x;
499 }
500 sh *= G4Log(tau/taul)/G4Log(taulim/taul);
501 }
502 sh *= 0.5;
503 return sh;
504 }
505
506
507 //....oooOO0OOooo........oooOO0OOooo........oo
508
509 G4double G4EmCorrections::ShellCorrection(cons
510 cons
511 cons
512 {
513 SetupKinematics(p, mat, ekin);
514 G4double term = 0.0;
515 //G4cout << "### G4EmCorrections::ShellCorre
516 // << " " << ekin/MeV << " MeV " <<
517 for (G4int i = 0; i<numberOfElements; ++i) {
518
519 G4double res = 0.0;
520 G4double res0 = 0.0;
521 const G4double Z = (*theElementVector)[i]-
522 const G4int iz = (*theElementVector)[i]->G
523 G4double Z2= (Z-0.3)*(Z-0.3);
524 G4double f = 1.0;
525 if(1 == iz) {
526 f = 0.5;
527 Z2 = 1.0;
528 }
529 G4double eta = ba2/Z2;
530 G4double tet = (11 < iz) ? sThetaK->Value(
531 res0 = f*KShell(tet,eta);
532 res += res0;
533 //G4cout << " Z= " << iz << " Shell 0" <<
534 // << " eta= " << eta << " resK= "
535
536 if(2 < iz) {
537 const G4double Zeff = (iz < 10) ? Z - ZD
538 Z2= Zeff*Zeff;
539 eta = ba2/Z2;
540 tet = sThetaL->Value(Z);
541 f = 0.125;
542 const G4int ntot = G4AtomicShells::GetNu
543 const G4int nmax = std::min(4, ntot);
544 G4double norm = 0.0;
545 G4double eshell = 0.0;
546 for(G4int j=1; j<nmax; ++j) {
547 G4int ne = G4AtomicShells::GetNumberOf
548 if(15 >= iz) {
549 tet = (3 > j) ? 0.25*Z2*(1.0 + 5*Z2*
550 0.25*Z2*(1.0 + Z2*alpha2/16.);
551 }
552 norm += ne;
553 eshell += tet*ne;
554 res0 = f*ne*LShell(tet,eta);
555 res += res0;
556 //G4cout << " Zeff= " << Zeff << " She
557 // << " tet= " << tet << " eta=
558 // << " resL= " << res0 << G4en
559 }
560 if(ntot > nmax) {
561 if (norm > 0.0) { norm = 1.0/norm; }
562 eshell *= norm;
563
564 static const G4double HM[53] = {
565 12.0, 12.0, 12.0, 12.0, 11.9, 11.7, 11.5,
566 10.0, 9.51, 8.97, 8.52, 8.03, 7.46, 6.95,
567 5.61, 5.39, 5.19, 5.01, 4.86, 4.72, 4.62,
568 4.32, 4.26, 4.20, 4.15, 4.1, 4.04, 4.00,
569 3.90, 3.89, 3.89, 3.88, 3.88, 3.88, 3.88,
570 3.90, 3.92, 3.93 };
571 static const G4double HN[31] = {
572 75.5, 61.9, 52.2, 45.1, 39.6, 35.4, 31.9,
573 24.5, 23.6, 22.7, 22.0, 21.4, 20.9, 20.5,
574 19.5, 19.3, 19.2, 19.1, 18.4, 18.8, 18.7,
575 18.2};
576
577 // Add M-shell
578 if(28 > iz) {
579 res += f*(iz - 10)*LShell(eshell,HM[
580 } else if(63 > iz) {
581 res += f*18*LShell(eshell,HM[iz-11]*
582 } else {
583 res += f*18*LShell(eshell,HM[52]*eta
584 }
585 // Add N-shell
586 if(32 < iz) {
587 if(60 > iz) {
588 res += f*(iz - 28)*LShell(eshell,H
589 } else if(63 > iz) {
590 res += 4*LShell(eshell,HN[iz-33]*e
591 } else {
592 res += 4*LShell(eshell,HN[30]*eta)
593 }
594 // Add O-P-shells
595 if(60 < iz) {
596 res += f*(iz - 60)*LShell(eshell,1
597 }
598 }
599 }
600 }
601 term += res*atomDensity[i]/Z;
602 }
603
604 term /= material->GetTotNbOfAtomsPerVolume()
605 //G4cout << "##Shell Correction=" << term <<
606 return term;
607 }
608
609 //....oooOO0OOooo........oooOO0OOooo........oo
610
611 G4double G4EmCorrections::DensityCorrection(co
612 co
613 co
614 {
615 SetupKinematics(p, mat, e);
616
617 G4double cden = material->GetIonisation()->
618 G4double mden = material->GetIonisation()->
619 G4double aden = material->GetIonisation()->
620 G4double x0den = material->GetIonisation()->
621 G4double x1den = material->GetIonisation()->
622
623 G4double dedx = 0.0;
624
625 // density correction
626 static const G4double twoln10 = 2.0*G4Log(10
627 G4double x = G4Log(bg2)/twoln10;
628 if ( x >= x0den ) {
629 dedx = twoln10*x - cden ;
630 if ( x < x1den ) dedx += aden*G4Exp(G4Log(
631 }
632
633 return dedx;
634 }
635
636 //....oooOO0OOooo........oooOO0OOooo........oo
637
638 G4double G4EmCorrections::BarkasCorrection(con
639 con
640 con
641 con
642 {
643 // . Z^3 Barkas effect in the stopping power
644 // J.C Ashley and R.H.Ritchie
645 // Physical review B Vol.5 No.7 1 April 19
646 // valid for kineticEnergy > 0.5 MeV
647
648 if (!isInitialized) { SetupKinematics(p, mat
649 G4double BarkasTerm = 0.0;
650
651 for (G4int i = 0; i<numberOfElements; ++i) {
652
653 const G4int iz = (*theElementVector)[i]->G
654 if(iz == 47) {
655 BarkasTerm += atomDensity[i]*0.006812*G4
656 } else if(iz >= 64) {
657 BarkasTerm += atomDensity[i]*0.002833*G4
658 } else {
659
660 const G4double Z = (*theElementVector)[i
661 const G4double X = ba2 / Z;
662 G4double b = 1.3;
663 if(1 == iz) { b = (material->GetName() =
664 else if(2 == iz) { b = 0.6; }
665 else if(10 >= iz) { b = 1.8; }
666 else if(17 >= iz) { b = 1.4; }
667 else if(18 == iz) { b = 1.8; }
668 else if(25 >= iz) { b = 1.4; }
669 else if(50 >= iz) { b = 1.35;}
670
671 const G4double W = b/std::sqrt(X);
672
673 G4double val = sBarkasCorr->Value(W, idx
674 if (W > sWmaxBarkas) { val *= (sWmaxBark
675 // G4cout << "i= " << i << " b= " <<
676 // << " Z= " << Z << " X= " << X << " va
677 BarkasTerm += val*atomDensity[i] / (std:
678 }
679 }
680
681 BarkasTerm *= 1.29*charge/material->GetTotNb
682
683 return BarkasTerm;
684 }
685
686 //....oooOO0OOooo........oooOO0OOooo........oo
687
688 G4double G4EmCorrections::BlochCorrection(cons
689 cons
690 cons
691 cons
692 {
693 if (!isInitialized) { SetupKinematics(p, mat
694
695 G4double y2 = q2/ba2;
696
697 G4double term = 1.0/(1.0 + y2);
698 G4double del;
699 G4double j = 1.0;
700 do {
701 j += 1.0;
702 del = 1.0/(j* (j*j + y2));
703 term += del;
704 // Loop checking, 03-Aug-2015, Vladimir Iv
705 } while (del > 0.01*term);
706
707 return -y2*term;
708 }
709
710 //....oooOO0OOooo........oooOO0OOooo........oo
711
712 G4double G4EmCorrections::MottCorrection(const
713 const
714 const
715 const
716 {
717 if (!isInitialized) { SetupKinematics(p, mat
718 return CLHEP::pi*CLHEP::fine_structure_const
719 }
720
721 //....oooOO0OOooo........oooOO0OOooo........oo
722
723 G4double
724 G4EmCorrections::EffectiveChargeCorrection(con
725 con
726 con
727 {
728 G4double factor = 1.0;
729 if(p->GetPDGCharge() <= 2.5*CLHEP::eplus ||
730
731 if(verbose > 1) {
732 G4cout << "EffectiveChargeCorrection: " <<
733 << " in " << mat->GetName()
734 << " ekin(MeV)= " << ekin/MeV << G4
735 }
736
737 if(p != curParticle || mat != curMaterial) {
738 curParticle = p;
739 curMaterial = mat;
740 curVector = nullptr;
741 currentZ = p->GetAtomicNumber();
742 if(verbose > 1) {
743 G4cout << "G4EmCorrections::EffectiveCha
744 << currentZ << " Aion= " << p->Ge
745 }
746 massFactor = CLHEP::proton_mass_c2/p->GetP
747 idx = -1;
748
749 for(G4int i=0; i<nIons; ++i) {
750 if(materialList[i] == mat && currentZ ==
751 idx = i;
752 break;
753 }
754 }
755 //G4cout << " idx= " << idx << " dz= " <<
756 if(idx >= 0) {
757 if(nullptr == ionList[idx]) { BuildCorre
758 curVector = stopData[idx];
759 } else {
760 return factor;
761 }
762 }
763 if(nullptr != curVector) {
764 factor = curVector->Value(ekin*massFactor)
765 if(verbose > 1) {
766 G4cout << "E= " << ekin << " factor= " <
767 << massFactor << G4endl;
768 }
769 }
770 return factor;
771 }
772
773 //....oooOO0OOooo........oooOO0OOooo........oo
774
775 void G4EmCorrections::AddStoppingData(const G4
776 const G4
777 G4Physic
778 {
779 G4int i = 0;
780 for(; i<nIons; ++i) {
781 if(Z == Zion[i] && A == Aion[i] && mname =
782 }
783 if(i == nIons) {
784 Zion.push_back(Z);
785 Aion.push_back(A);
786 materialName.push_back(mname);
787 materialList.push_back(nullptr);
788 ionList.push_back(nullptr);
789 stopData.push_back(dVector);
790 nIons++;
791 if(verbose > 1) {
792 G4cout << "AddStoppingData Z= " << Z <<
793 << " idx= " << i << G4endl;
794 }
795 }
796 }
797
798 //....oooOO0OOooo........oooOO0OOooo........oo
799
800 void G4EmCorrections::BuildCorrectionVector()
801 {
802 if(nullptr == ionLEModel || nullptr == ionHE
803 return;
804 }
805
806 const G4ParticleDefinition* ion = curParticl
807 const G4ParticleDefinition* gion = G4Generic
808 G4int Z = Zion[idx];
809 G4double A = Aion[idx];
810 G4PhysicsVector* v = stopData[idx];
811
812 if(verbose > 1) {
813 G4cout << "BuildCorrectionVector: Stopping
814 << curParticle->GetParticleName() <
815 << materialName[idx] << " Ion Z= "
816 << " massFactor= " << massFactor <<
817 G4cout << " Nbins=" << nbinCorr << " Em
818 << " Emax(MeV)=" << eCorrMax << " ion: "
819 << ion->GetParticleName() << G4endl
820 }
821
822 auto vv = new G4PhysicsLogVector(eCorrMin,eC
823 const G4double eth0 = v->Energy(0);
824 const G4double escal = eth/massFactor;
825 G4double qe =
826 effCharge.EffectiveChargeSquareRatio(curPa
827 const G4double dedxt =
828 ionLEModel->ComputeDEDXPerVolume(curMateri
829 const G4double dedx1t =
830 ionHEModel->ComputeDEDXPerVolume(curMateri
831 + ComputeIonCorrections(curParticle, curMa
832 const G4double rest = escal*(dedxt - dedx1t)
833 if(verbose > 1) {
834 G4cout << "Escal(MeV)= " << escal << " qe=
835 << " dedxt= " << dedxt << " dedx1t=
836 }
837 for(G4int i=0; i<=nbinCorr; ++i) {
838 // energy in the local table (GenericIon)
839 G4double e = vv->Energy(i);
840 // energy of the real ion
841 G4double eion = e/massFactor;
842 // energy in the imput stopping data vecto
843 G4double e1 = eion/A;
844 G4double dedx = (e1 < eth0)
845 ? v->Value(eth0)*std::sqrt(e1/eth0) : v-
846 qe = effCharge.EffectiveChargeSquareRatio(
847 G4double dedx1 = (e <= eth)
848 ? ionLEModel->ComputeDEDXPerVolume(curMa
849 : ionHEModel->ComputeDEDXPerVolume(curMa
850 ComputeIonCorrections(curParticle, curMa
851 vv->PutValue(i, dedx/dedx1);
852 if(verbose > 1) {
853 G4cout << "E(MeV)=" << e/CLHEP::MeV << "
854 << " e1=" << e1 << " dedxRatio= " << de
855 << " dedx=" << dedx << " dedx1="
856 << " qe=" << qe << " rest/eion="
857 }
858 }
859 delete v;
860 ionList[idx] = ion;
861 stopData[idx] = vv;
862 if(verbose > 1) { G4cout << "End data set "
863 }
864
865 //....oooOO0OOooo........oooOO0OOooo........oo
866
867 void G4EmCorrections::InitialiseForNewRun()
868 {
869 G4ProductionCutsTable* tb = G4ProductionCuts
870 ncouples = tb->GetTableSize();
871 if(currmat.size() != ncouples) {
872 currmat.resize(ncouples);
873 for(auto it = thcorr.begin(); it != thcorr
874 (it->second).clear();
875 }
876 thcorr.clear();
877 for(std::size_t i=0; i<ncouples; ++i) {
878 currmat[i] = tb->GetMaterialCutsCouple((
879 G4String nam = currmat[i]->GetName();
880 for(G4int j=0; j<nIons; ++j) {
881 if(nam == materialName[j]) { materialL
882 }
883 }
884 }
885 }
886
887 //....oooOO0OOooo........oooOO0OOooo........oo
888
889 void G4EmCorrections::Initialise()
890 {
891 // Z^3 Barkas effect in the stopping power o
892 // J.C Ashley and R.H.Ritchie
893 // Physical review B Vol.5 No.7 1 April 1972
894 // G.S. Khandelwal Nucl. Phys. A116(1968)97
895 // "Shell corrections for K- and L- electron
896
897 G4int i, j;
898 static const G4double fTable[47][2] = {
899 { 0.02, 21.5},
900 { 0.03, 20.0},
901 { 0.04, 18.0},
902 { 0.05, 15.6},
903 { 0.06, 15.0},
904 { 0.07, 14.0},
905 { 0.08, 13.5},
906 { 0.09, 13.},
907 { 0.1, 12.2},
908 { 0.2, 9.25},
909 { 0.3, 7.0},
910 { 0.4, 6.0},
911 { 0.5, 4.5},
912 { 0.6, 3.5},
913 { 0.7, 3.0},
914 { 0.8, 2.5},
915 { 0.9, 2.0},
916 { 1.0, 1.7},
917 { 1.2, 1.2},
918 { 1.3, 1.0},
919 { 1.4, 0.86},
920 { 1.5, 0.7},
921 { 1.6, 0.61},
922 { 1.7, 0.52},
923 { 1.8, 0.5},
924 { 1.9, 0.43},
925 { 2.0, 0.42},
926 { 2.1, 0.3},
927 { 2.4, 0.2},
928 { 3.0, 0.13},
929 { 3.08, 0.1},
930 { 3.1, 0.09},
931 { 3.3, 0.08},
932 { 3.5, 0.07},
933 { 3.8, 0.06},
934 { 4.0, 0.051},
935 { 4.1, 0.04},
936 { 4.8, 0.03},
937 { 5.0, 0.024},
938 { 5.1, 0.02},
939 { 6.0, 0.013},
940 { 6.5, 0.01},
941 { 7.0, 0.009},
942 { 7.1, 0.008},
943 { 8.0, 0.006},
944 { 9.0, 0.0032},
945 { 10.0, 0.0025} };
946
947 sBarkasCorr = new G4PhysicsFreeVector(47, fa
948 for(i=0; i<47; ++i) { sBarkasCorr->PutValues
949
950 const G4double SK[20] = {1.9477, 1.9232, 1.8
951 1.7754, 1.7396, 1.7
952 1.6461, 1.6189, 1.5
953 1.5467, 1.5254, 1.5
954 const G4double TK[20] = {2.5222, 2.5125, 2.5
955 2.4388, 2.4163, 2.4
956 2.3466, 2.3229, 2.2
957 2.2515, 2.2277, 2.2
958
959 const G4double SL[26] = {15.3343, 13.9389, 1
960 10.3424, 10.0371,
961 8.4114, 8.0683,
962 7.2506, 7.0327,
963 6.4969, 6.3498,
964 const G4double TL[26] = {35.0669, 33.4344, 3
965 28.6128, 28.1449, 2
966 25.4058, 24.7587, 2
967 23.0771, 22.5880, 2
968 21.2872, 20.9006, 2
969
970 const G4double bk1[29][11] = {
971 {0.005, 1.34782E-8, 1.46132E-8, 1.72179E-8,
972 {0.007, 6.87555E-8, 7.44373E-8, 8.74397E-8,
973 {0.01, 3.78413E-7, 4.08831E-7, 4.78154E-7, 5
974 {0.015, 2.53200E-6, 2.72664E-6, 3.16738E-6,
975 {0.02, 9.43891E-6, 1.01339E-5, 1.16984E-5, 1
976 {0.03, 5.67088E-5, 6.05576E-5, 6.91311E-5, 7
977 {0.04, 1.91576E-4, 2.03626E-4, 2.30230E-4, 2
978 {0.05, 4.74226E-4, 5.02006E-4, 5.62872E-4, 6
979 {0.06, 9.67285E-4, 1.02029E-3, 1.13566E-3, 1
980 {0.08, 2.81537E-3, 2.95200E-3, 3.24599E-3, 3
981 {0.1, 6.14216E-3, 6.40864E-3, 6.97750E-3, 7.
982 {0.15, 2.23096E-2, 2.30790E-2, 2.46978E-2, 2
983 {0.2, 5.04022E-2, 5.18492E-2, 5.48682E-2, 5.
984 {0.3, 1.38018E-1, 1.41026E-1, 1.47244E-1, 1.
985 {0.4, 2.56001E-1, 2.60576E-1, 2.69992E-1, 2.
986 {0.5, 3.92181E-1, 3.98213E-1, 4.10597E-1, 4.
987 {0.6, 5.37913E-1, 5.45268E-1, 5.60350E-1, 5.
988 {0.7, 6.87470E-1, 6.96021E-1, 7.13543E-1, 7.
989 {0.8, 8.37159E-1, 8.46790E-1, 8.66525E-1, 8.
990 {1.0, 1.12850, 1.14002, 1.16362, 1.18799, 1.
991 {1.2, 1.40232, 1.41545, 1.44232, 1.47007, 1.
992 {1.4, 1.65584, 1.67034, 1.70004, 1.73072, 1.
993 {1.5, 1.77496, 1.79009, 1.82107, 1.85307, 1.
994 {1.7, 1.99863, 2.01490, 2.04822, 2.08265, 2.
995 {2.0, 2.29325, 2.31100, 2.34738, 2.38501, 2.
996 {3.0, 3.08887, 3.11036, 3.15445, 3.20013, 3.
997 {5.0, 4.07599, 4.10219, 4.15606, 4.21199, 4.
998 {7.0, 4.69647, 4.72577, 4.78608, 4.84877, 4.
999 {10.0, 5.32590, 5.35848, 5.42560, 5.49547, 5
1000 };
1001
1002 const G4double bk2[29][11] = {
1003 {0.005, 5.98040E-8, 7.25636E-8, 8.00602E-8,
1004 {0.007, 2.95345E-7, 3.56497E-7, 3.92247E-7,
1005 {0.01, 1.55232E-6, 1.86011E-6, 2.03881E-6,
1006 {0.015, 9.67802E-6, 1.14707E-5, 1.25008E-5,
1007 {0.02, 3.38425E-5, 3.97259E-5, 4.30763E-5,
1008 {0.03, 1.81920E-4, 2.10106E-4, 2.25918E-4,
1009 {0.04, 5.59802E-4, 6.38103E-4, 6.81511E-4,
1010 {0.05, 1.28002E-3, 1.44336E-3, 1.53305E-3,
1011 {0.06, 2.42872E-3, 2.72510E-3, 2.88111E-3,
1012 {0.08, 6.35222E-3, 6.99730E-3, 7.34446E-3,
1013 {0.1, 1.26929E-2, 1.38803E-2, 1.44371E-2, 1
1014 {0.15, 3.96872E-2, 4.24699E-2, 4.39340E-2,
1015 {0.2, 8.15290E-2, 8.62830E-2, 8.87650E-2, 9
1016 {0.3, 1.99528E-1, 2.08471E-1, 2.13103E-1, 2
1017 {0.4, 3.47383E-1, 3.60369E-1, 3.67073E-1, 3
1018 {0.5, 5.11214E-1, 5.27935E-1, 5.36554E-1, 5
1019 {0.6, 6.82122E-1, 7.02260E-1, 7.12631E-1, 7
1020 {0.7, 8.54544E-1, 8.77814E-1, 8.89791E-1, 9
1021 {0.8, 1.02508, 1.05121, 1.06466, 1.07838, 1
1022 {1.0, 1.35307, 1.38429, 1.40036, 1.41676, 1
1023 {1.2, 1.65823, 1.69385, 1.71220, 1.73092, 1
1024 {1.4, 1.93902, 1.97852, 1.99887, 2.01964, 2
1025 {1.5, 2.07055, 2.11182, 2.13309, 2.15480, 2
1026 {1.7, 2.31700, 2.36154, 2.38451, 2.40798, 2
1027 {2.0, 2.64162, 2.69053, 2.71576, 2.74154, 2
1028 {3.0, 3.51376, 3.57394, 3.60503, 3.63684, 3
1029 {5.0, 4.59935, 4.67433, 4.71316, 4.75293, 4
1030 {7.0, 5.28542, 5.37040, 5.41445, 5.45962, 5
1031 {10.0, 5.98474, 6.08046, 6.13015, 6.18112,
1032 };
1033
1034 const G4double bls1[28][10] = {
1035 {0.005, 2.4111E-4, 2.5612E-4, 2.7202E-4, 3.
1036 {0.007, 6.3947E-4, 6.7058E-4, 7.0295E-4, 7.
1037 {0.01, 1.5469E-3, 1.6036E-3, 1.6622E-3, 1.7
1038 {0.015, 3.7221E-3, 3.8404E-3, 3.9650E-3, 4.
1039 {0.02, 6.9449E-3, 7.1910E-3, 7.4535E-3, 8.0
1040 {0.03, 1.7411E-2, 1.8180E-2, 1.8997E-2, 2.0
1041 {0.04, 3.8474E-2, 4.0056E-2, 4.1718E-2, 4.5
1042 {0.05, 6.7371E-2, 6.9928E-2, 7.2596E-2, 7.8
1043 {0.06, 1.0418E-1, 1.0778E-1, 1.1152E-1, 1.1
1044 {0.08, 1.9647E-1, 2.0217E-1, 2.0805E-1, 2.2
1045 {0.1, 3.0594E-1, 3.1361E-1, 3.2150E-1, 3.37
1046 {0.15, 6.1411E-1, 6.2597E-1, 6.3811E-1, 6.6
1047 {0.2, 9.3100E-1, 9.5614E-1, 9.7162E-1, 1.00
1048 {0.3, 1.5172, 1.5372, 1.5576, 1.5998, 1.643
1049 {0.4, 2.0173, 2.0408, 2.0647, 2.1142, 2.165
1050 {0.5, 2.3932, 2.4193, 2.4460, 2.5011, 2.558
1051 {0.6, 2.7091, 2.7374, 2.7663, 2.8260, 2.888
1052 {0.7, 2.9742, 3.0044, 3.0352, 3.0988, 3.165
1053 {0.8, 3.2222, 3.2539, 3.2863, 3.3532, 3.423
1054 {1.0, 3.6690, 3.7033, 3.7384, 3.8108, 3.886
1055 {1.2, 3.9819, 4.0183, 4.0555, 4.1324, 4.213
1056 {1.4, 4.2745, 4.3127, 4.3517, 4.4324, 4.517
1057 {1.5, 4.4047, 4.4436, 4.4834, 4.5658, 4.652
1058 {1.7, 4.6383, 4.6787, 4.7200, 4.8054, 4.894
1059 {2.0, 4.9369, 4.9791, 5.0223, 5.1116, 5.205
1060 {3.0, 5.6514, 5.6981, 5.7459, 5.8450, 5.948
1061 {5.0, 6.4665, 6.5189, 6.5724, 6.6835, 6.800
1062 {7.0, 6.8634, 6.9194, 6.9767, 7.0957, 7.220
1063 };
1064
1065 const G4double bls2[28][10] = {
1066 {0.005, 5.4561E-4, 6.0905E-4, 6.7863E-4, 7.
1067 {0.007, 1.2068E-3, 1.3170E-3, 1.4377E-3, 1.
1068 {0.01, 2.6832E-3, 2.9017E-3, 3.1534E-3, 3.4
1069 {0.015, 6.2775E-3, 6.9077E-3, 7.6525E-3, 8.
1070 {0.02, 1.2561E-2, 1.3943E-2, 1.5553E-2, 1.7
1071 {0.03, 3.3750E-2, 3.7470E-2, 4.1528E-2, 4.6
1072 {0.04, 6.9619E-2, 7.6098E-2, 8.3249E-2, 9.1
1073 {0.05, 1.1522E-1, 1.2470E-1, 1.3504E-1, 1.4
1074 {0.06, 1.6931E-1, 1.8179E-1, 1.9530E-1, 2.0
1075 {0.08, 2.9540E-1, 3.1361E-1, 3.3312E-1, 3.5
1076 {0.1, 4.3613E-1, 4.5956E-1, 4.852E-1, 5.1
1077 {0.15, 8.1014E-1, 8.4453E-1, 8.8093E-1, 9.1
1078 {0.2, 1.1888, 1.2319, 1.2774, 1.3255, 1.350
1079 {0.3, 1.8422, 1.8983, 1.9575, 2.0201, 2.052
1080 {0.4, 2.3984, 2.4642, 2.5336, 2.6070, 2.645
1081 {0.5, 2.8181, 2.8915, 2.9690, 3.0509, 3.093
1082 {0.6, 3.1698, 3.2494, 3.3336, 3.4226, 3.469
1083 {0.7, 3.4652, 3.5502, 3.6400, 3.7351, 3.784
1084 {0.8, 3.7392, 3.8289, 3.9236, 4.0239, 4.076
1085 {1.0, 4.2295, 4.3269, 4.4299, 4.5391, 4.596
1086 {1.2, 4.5777, 4.6814, 4.7912, 4.9076, 4.968
1087 {1.4, 4.9001, 5.0092, 5.1247, 5.2473, 5.311
1088 {1.5, 5.0434, 5.1550, 5.2731, 5.3984, 5.464
1089 {1.7, 5.3011, 5.4170, 5.5398, 5.6701, 5.737
1090 {2.0, 5.6308, 5.7523, 5.8811, 6.0174, 6.089
1091 {3.0, 6.4224, 6.5580, 6.7019, 6.8549, 6.935
1092 {5.0, 7.3338, 7.4872, 7.6500, 7.8235, 7.914
1093 {7.0, 7.7938, 7.9588, 8.1342, 8.3211, 8.419
1094 };
1095
1096 const G4double bls3[28][9] = {
1097 {0.005, 1.2895E-3, 1.3670E-3, 1.4524E-3, 1.
1098 {0.007, 2.6467E-3, 2.8242E-3, 3.0238E-3, 3.
1099 {0.01, 6.1472E-3, 6.6086E-3, 7.1246E-3, 8.3
1100 {0.015, 1.63316E-2, 1.7572E-2, 1.8932E-2, 2
1101 {0.02, 3.2634E-2, 3.4900E-2, 3.7348E-2, 4.2
1102 {0.03, 7.9907E-2, 8.4544E-2, 8.9486E-2, 1.0
1103 {0.04, 1.4523E-1, 1.5237E-1, 1.5985E-1, 1.7
1104 {0.05, 2.2082E-1, 2.3036E-1, 2.4038E-1, 2.6
1105 {0.06, 3.0423E-1, 3.1611E-1, 3.2854E-1, 3.5
1106 {0.08, 4.8536E-1, 5.0156E-1, 5.1846E-1, 5.5
1107 {0.1, 6.7586E-1, 6.9596E-1, 7.1688E-1, 7.61
1108 {0.15, 1.1544, 1.1828, 1.2122, 1.2746, 1.34
1109 {0.2, 1.6167, 1.6517, 1.6880, 1.7650, 1.848
1110 {0.3, 2.3979, 2.4432, 2.4902, 2.5899, 2.698
1111 {0.4, 3.0502, 3.1034, 3.1586, 3.2758, 3.403
1112 {0.5, 3.5466, 3.6062, 3.6681, 3.7994, 3.942
1113 {0.6, 3.9620, 4.0270, 4.0945, 4.2378, 4.393
1114 {0.7, 4.3020, 4.3715, 4.4438, 4.5974, 4.764
1115 {0.8, 4.6336, 4.7072, 4.7837, 4.9463, 5.123
1116 {1.0, 5.2041, 5.2845, 5.3682, 5.5462, 5.740
1117 {1.2, 5.6182, 5.7044, 5.7940, 5.9848, 6.192
1118 {1.4, 5.9967, 6.0876, 6.1823, 6.3839, 6.603
1119 {1.5, 6.1652, 6.2583, 6.3553, 6.5618, 6.787
1120 {1.7, 6.4686, 6.5657, 6.6668, 6.8823, 7.117
1121 {2.0, 6.8577, 6.9600, 7.0666, 7.2937, 7.541
1122 {3.0, 7.7981, 7.9134, 8.0336, 8.2901, 8.570
1123 {5.0, 8.8978, 9.0297, 9.1673, 9.4612, 9.783
1124 {7.0, 9.4819, 9.6248, 9.7739, 10.0926, 10.4
1125 };
1126
1127 const G4double bll1[28][10] = {
1128 {0.005, 3.6324E-5, 4.0609E-5, 4.5430E-5, 5.
1129 {0.007, 1.8110E-4, 2.0001E-4, 2.2099E-4, 2.
1130 {0.01, 8.6524E-4, 9.4223E-4, 1.0262E-3, 1.2
1131 {0.015, 4.2293E-3, 4.5314E-3, 4.8551E-3, 5.
1132 {0.02, 1.1485E-2, 1.2172E-2, 1.2900E-2, 1.4
1133 {0.03, 3.9471E-2, 4.1270E-2, 4.3149E-2, 4.7
1134 {0.04, 8.4454E-2, 8.7599E-2, 9.0860E-2, 9.7
1135 {0.05, 1.4339E-1, 1.4795E-1, 1.5266E-1, 1.6
1136 {0.06, 2.1304E-1, 2.1899E-1, 2.2512E-1, 2.3
1137 {0.08, 3.7382E-1, 3.8241E-1, 3.9122E-1, 4.0
1138 {0.1, 5.5056E-1, 5.6151E-1, 5.7273E-1, 5.96
1139 {0.15, 1.0066, 1.0224, 1.0386, 1.0721, 1.10
1140 {0.2, 1.4376, 1.4572, 1.4773, 1.5188, 1.562
1141 {0.3, 2.1712, 2.1964, 2.2223, 2.2758, 2.332
1142 {0.4, 2.7500, 2.7793, 2.8094, 2.8719, 2.937
1143 {0.5, 3.2033, 3.2359, 3.2693, 3.3389, 3.412
1144 {0.6, 3.6038, 3.6391, 3.6753, 3.7506, 3.830
1145 {0.7, 3.9106, 3.9482, 3.9867, 4.0670, 4.152
1146 {0.8, 4.1790, 4.2185, 4.2590, 4.3437, 4.433
1147 {1.0, 4.6344, 4.6772, 4.7212, 4.8131, 4.910
1148 {1.2, 4.9787, 5.0242, 5.0711, 5.1689, 5.272
1149 {1.4, 5.2688, 5.3166, 5.3658, 5.4687, 5.578
1150 {1.5, 5.3966, 5.4454, 5.4957, 5.6009, 5.712
1151 {1.7, 5.6255, 5.6762, 5.7284, 5.8377, 5.954
1152 {2.0, 5.9170, 5.9701, 6.0248, 6.1395, 6.261
1153 {3.0, 6.6210, 6.6801, 6.7411, 6.8692, 7.006
1154 {5.0, 7.4620, 7.5288, 7.5977, 7.7428, 7.898
1155 {7.0, 7.7362, 7.8079, 7.8821, 8.0383, 8.206
1156 };
1157
1158 const G4double bll2[28][10] = {
1159 {0.005, 1.8339E-4, 2.3330E-4, 2.9738E-4, 3.
1160 {0.007, 7.5042E-4, 9.2355E-4, 1.1375E-3, 1.
1161 {0.01, 2.8829E-3, 3.4275E-3, 4.0758E-3, 4.8
1162 {0.015, 1.1087E-2, 1.2716E-2, 1.4581E-2, 1.
1163 {0.02, 2.5786E-2, 2.8922E-2, 3.2435E-2, 3.6
1164 {0.03, 7.3461E-2, 8.0264E-2, 8.7705E-2, 9.5
1165 {0.04, 1.4094E-1, 1.5172E-1, 1.6336E-1, 1.7
1166 {0.05, 2.2289E-1, 2.3762E-1, 2.5344E-1, 2.7
1167 {0.06, 3.1503E-1, 3.3361E-1, 3.5346E-1, 3.7
1168 {0.08, 5.1793E-1, 5.4368E-1, 5.7109E-1, 6.0
1169 {0.1, 7.3258E-1, 7.6481E-1, 7.9907E-1, 8.35
1170 {0.15, 1.2677, 1.3137, 1.3626, 1.4147, 1.44
1171 {0.2, 1.7615, 1.8188, 1.8797, 1.9446, 1.978
1172 {0.3, 2.5909, 2.6658, 2.7457, 2.8312, 2.876
1173 {0.4, 3.2417, 3.3302, 3.4249, 3.5266, 3.580
1174 {0.5, 3.7527, 3.8523, 3.9591, 4.0741, 4.135
1175 {0.6, 4.2013, 4.3103, 4.4274, 4.5537, 4.620
1176 {0.7, 4.5493, 4.6664, 4.7925, 4.9286, 5.000
1177 {0.8, 4.8537, 4.9780, 5.1119, 5.2568, 5.333
1178 {1.0, 5.3701, 5.5066, 5.6540, 5.8138, 5.898
1179 {1.2, 5.7648, 5.9114, 6.0701, 6.2424, 6.334
1180 {1.4, 6.0976, 6.2530, 6.4214, 6.6044, 6.702
1181 {1.5, 6.2447, 6.4041, 6.5768, 6.7647, 6.865
1182 {1.7, 6.5087, 6.6752, 6.8558, 7.0526, 7.157
1183 {2.0, 6.8458, 7.0218, 7.2129, 7.4213, 7.532
1184 {3.0, 7.6647, 7.8644, 8.0819, 8.3189, 8.447
1185 {5.0, 8.6515, 8.8816, 9.1330, 9.4090, 9.557
1186 {7.0, 9.0221, 9.2724, 9.5464, 9.8477, 10.00
1187 };
1188
1189 const G4double bll3[28][9] = {
1190 {0.005, 1.3190E-3, 1.4961E-3, 1.6974E-3, 2.
1191 {0.007, 4.0158E-3, 4.4623E-3, 4.9592E-3, 6.
1192 {0.01, 1.1509E-2, 1.2548E-2, 1.3681E-2, 1.6
1193 {0.015, 3.3070E-2, 3.5408E-2, 3.7914E-2, 4.
1194 {0.02, 6.4555E-2, 6.8394E-2, 7.2472E-2, 8.1
1195 {0.03, 1.5030E-1, 1.5101E-1, 1.5844E-1, 1.7
1196 {0.04, 2.5743E-1, 2.6774E-1, 2.7855E-1, 3.0
1197 {0.05, 3.7846E-1, 3.9195E-1, 4.0607E-1, 4.3
1198 {0.06, 5.0839E-1, 5.2497E-1, 5.4230E-1, 5.7
1199 {0.08, 7.8230E-1, 8.0474E-1, 8.2818E-1, 8.7
1200 {0.1, 1.0621, 1.0900, 1.1192, 1.1816, 1.250
1201 {0.15, 1.7389, 1.7790, 1.8210, 1.9112, 2.01
1202 {0.2, 2.3516, 2.4024, 2.4556, 2.5701, 2.697
1203 {0.3, 3.3741, 3.4427, 3.5148, 3.6706, 3.844
1204 {0.4, 4.1788, 4.2620, 4.3496, 4.5398, 4.753
1205 {0.5, 4.8180, 4.9137, 5.0146, 5.2341, 5.481
1206 {0.6, 5.3765, 5.4830, 5.5954, 5.8406, 6.117
1207 {0.7, 5.8208, 5.9369, 6.0596, 6.3275, 6.630
1208 {0.8, 6.2109, 6.3355, 6.4674, 6.7558, 7.082
1209 {1.0, 6.8747, 7.0142, 7.1621, 7.4861, 7.854
1210 {1.2, 7.3933, 7.5454, 7.7068, 8.0612, 8.465
1211 {1.4, 7.8331, 7.9967, 8.1694, 8.5502, 8.985
1212 {1.5, 8.0286, 8.1967, 8.3753, 8.7681, 9.218
1213 {1.7, 8.3813, 8.5585, 8.7469, 9.1618, 9.636
1214 {2.0, 8.8352, 9.0245, 9.2260, 9.6701, 10.17
1215 {3.0, 9.9511, 10.1714, 10.4062, 10.9254, 11
1216 {5.0, 11.3211, 11.5818, 11.8601, 12.4771, 1
1217 {7.0, 11.9480, 12.2357, 12.5432, 13.2260, 1
1218 };
1219
1220 G4double b, bs;
1221 for(i=0; i<nEtaK; ++i) {
1222
1223 G4double et = Eta[i];
1224 G4double loget = G4Log(et);
1225
1226 for(j=0; j<nK; ++j) {
1227
1228 if(j < 10) { b = bk2[i][10-j]; }
1229 else { b = bk1[i][20-j]; }
1230
1231 CK[j][i] = SK[j]*loget + TK[j] - b;
1232
1233 if(i == nEtaK-1) {
1234 ZK[j] = et*(et*et*CK[j][i] - et*UK[j]
1235 //G4cout << "i= " << i << " j= " << j
1236 // << " CK[j][i]= " << CK[j][i
1237 // << " ZK[j]= " << ZK[j] << "
1238 }
1239 }
1240 //G4cout << G4endl;
1241 if(i < nEtaL) {
1242 //G4cout << " LShell:" <<G4endl;
1243 for(j=0; j<nL; ++j) {
1244
1245 if(j < 8) {
1246 bs = bls3[i][8-j];
1247 b = bll3[i][8-j];
1248 } else if(j < 17) {
1249 bs = bls2[i][17-j];
1250 b = bll2[i][17-j];
1251 } else {
1252 bs = bls1[i][26-j];
1253 b = bll1[i][26-j];
1254 }
1255 G4double c = SL[j]*loget + TL[j];
1256 CL[j][i] = c - bs - 3.0*b;
1257 if(i == nEtaL-1) {
1258 VL[j] = et*(et*CL[j][i] - UL[j]);
1259 //G4cout << "i= " << i << " j= " <<
1260 // << " CL[j][i]= " << CL[
1261 // << " VL[j]= " << VL[j] <
1262 // << " et= " << et << G4en
1263 //" UL= " << UL[j] << " TL= " << TL
1264 }
1265 }
1266 //G4cout << G4endl;
1267 }
1268 }
1269
1270 const G4double xzk[34] = { 11.7711,
1271 13.3669, 15.5762, 17.1715, 18.7667, 20.85
1272 34.3457, 37.4119, 40.3555, 42.3177, 44.77
1273 60.9586, 63.6567, 66.5998, 68.807, 71.872
1274 93.4549, 96.2753
1275 const G4double yzk[34] = { 0.70663,
1276 0.72033, 0.73651, 0.74647, 0.75518, 0.763
1277 0.80611, 0.8123, 0.8185, 0.82097, 0.82467
1278 0.84565, 0.84936, 0.85181, 0.85303, 0.855
1279 0.86891, 0.87011
1280
1281 const G4double xzl[36] = { 15.5102,
1282 16.7347, 17.9592, 19.551, 21.0204, 22.612
1283 34.4898, 36.2041, 38.4082, 40.3674, 42.57
1284 59.3469, 61.9184, 64.6122, 67.4286, 71.46
1285 91.551, 94.2449,
1286 const G4double yzl[36] = { 0.29875,
1287 0.31746, 0.33368, 0.35239, 0.36985, 0.387
1288 0.4896, 0.50083, 0.51331, 0.52328, 0.5307
1289 0.58191, 0.5869, 0.59189, 0.60062, 0.6068
1290 0.6368, 0.64054
1291
1292 sThetaK = new G4PhysicsFreeVector(34, false
1293 for(i=0; i<34; ++i) { sThetaK->PutValues(i,
1294 sThetaL = new G4PhysicsFreeVector(36, false
1295 for(i=0; i<36; ++i) { sThetaL->PutValues(i,
1296 }
1297
1298 //....oooOO0OOooo........oooOO0OOooo........o
1299
1300