| 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 file
30 //
31 //
32 // File name: G4EmCalculator
33 //
34 // Author: Vladimir Ivanchenko
35 //
36 // Creation date: 28.06.2004
37 //
38 //
39 // Class Description: V.Ivanchenko & M.Novak
40 //
41 // -------------------------------------------
42 //
43 //....oooOO0OOooo........oooOO0OOooo........oo
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 #include "G4EmCalculator.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "G4LossTableManager.hh"
49 #include "G4EmParameters.hh"
50 #include "G4NistManager.hh"
51 #include "G4DynamicParticle.hh"
52 #include "G4VEmProcess.hh"
53 #include "G4VEnergyLossProcess.hh"
54 #include "G4VMultipleScattering.hh"
55 #include "G4Material.hh"
56 #include "G4MaterialCutsCouple.hh"
57 #include "G4ParticleDefinition.hh"
58 #include "G4ParticleTable.hh"
59 #include "G4IonTable.hh"
60 #include "G4PhysicsTable.hh"
61 #include "G4ProductionCutsTable.hh"
62 #include "G4ProcessManager.hh"
63 #include "G4ionEffectiveCharge.hh"
64 #include "G4RegionStore.hh"
65 #include "G4Element.hh"
66 #include "G4EmCorrections.hh"
67 #include "G4GenericIon.hh"
68 #include "G4ProcessVector.hh"
69 #include "G4Gamma.hh"
70 #include "G4Electron.hh"
71 #include "G4Positron.hh"
72 #include "G4EmUtility.hh"
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 G4EmCalculator::G4EmCalculator()
77 {
78 manager = G4LossTableManager::Instance();
79 nist = G4NistManager::Instance();
80 theParameters = G4EmParameters::Instance();
81 corr = manager->EmCorrections();
82 cutenergy[0] = cutenergy[1] = cutenergy[2] =
83 theGenericIon = G4GenericIon::GenericIon();
84 ionEffCharge = new G4ionEffectiveCharge();
85 dynParticle = new G4DynamicParticle();
86 ionTable = G4ParticleTable::GetParticle
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90
91 G4EmCalculator::~G4EmCalculator()
92 {
93 delete ionEffCharge;
94 delete dynParticle;
95 for (G4int i=0; i<nLocalMaterials; ++i) {
96 delete localCouples[i];
97 }
98 }
99
100 //....oooOO0OOooo........oooOO0OOooo........oo
101
102 G4double G4EmCalculator::GetDEDX(G4double kinE
103 const G4Parti
104 const G4Mater
105 const G4Regio
106 {
107 G4double res = 0.0;
108 const G4MaterialCutsCouple* couple = FindCou
109 if(nullptr != couple && UpdateParticle(p, ki
110 res = manager->GetDEDX(p, kinEnergy, coupl
111
112 if(isIon) {
113 if(FindEmModel(p, currentProcessName, ki
114 G4double length = CLHEP::nm;
115 G4double eloss = res*length;
116 //G4cout << "### GetDEDX: E= " << kinE
117 // << " de= " << eloss << G4endl
118 dynParticle->SetKineticEnergy(kinEnerg
119 currentModel->GetChargeSquareRatio(p,
120 currentModel->CorrectionsAlongStep(cou
121 res = eloss/length;
122 //G4cout << " de1= " << eloss <<
123 // << " " << p->GetParticleName(
124 }
125 }
126
127 if(verbose>0) {
128 G4cout << "G4EmCalculator::GetDEDX: E(Me
129 << " DEDX(MeV/mm)= " << res*mm/Me
130 << " DEDX(MeV*cm^2/g)= " << res*g
131 << " " << p->GetParticleName()
132 << " in " << mat->GetName()
133 << " isIon= " << isIon
134 << G4endl;
135 }
136 }
137 return res;
138 }
139
140 //....oooOO0OOooo........oooOO0OOooo........oo
141
142 G4double G4EmCalculator::GetRangeFromRestricte
143
144
145
146 {
147 G4double res = 0.0;
148 const G4MaterialCutsCouple* couple = FindCou
149 if(couple && UpdateParticle(p, kinEnergy)) {
150 res = manager->GetRangeFromRestricteDEDX(p
151 if(verbose>1) {
152 G4cout << " G4EmCalculator::GetRangeFrom
153 << kinEnergy/MeV
154 << " range(mm)= " << res/mm
155 << " " << p->GetParticleName()
156 << " in " << mat->GetName()
157 << G4endl;
158 }
159 }
160 return res;
161 }
162
163 //....oooOO0OOooo........oooOO0OOooo........oo
164
165 G4double G4EmCalculator::GetCSDARange(G4double
166 const G4
167 const G4
168 const G4
169 {
170 G4double res = 0.0;
171 if(!theParameters->BuildCSDARange()) {
172 G4ExceptionDescription ed;
173 ed << "G4EmCalculator::GetCSDARange: CSDA
174 << " use UI command: /process/eLoss/CSD
175 G4Exception("G4EmCalculator::GetCSDARange"
176 JustWarning, ed);
177 return res;
178 }
179
180 const G4MaterialCutsCouple* couple = FindCou
181 if(nullptr != couple && UpdateParticle(p, ki
182 res = manager->GetCSDARange(p, kinEnergy,
183 if(verbose>1) {
184 G4cout << " G4EmCalculator::GetCSDARange
185 << " range(mm)= " << res/mm
186 << " " << p->GetParticleName()
187 << " in " << mat->GetName()
188 << G4endl;
189 }
190 }
191 return res;
192 }
193
194 //....oooOO0OOooo........oooOO0OOooo........oo
195
196 G4double G4EmCalculator::GetRange(G4double kin
197 const G4Part
198 const G4Mate
199 const G4Regi
200 {
201 G4double res = 0.0;
202 if(theParameters->BuildCSDARange()) {
203 res = GetCSDARange(kinEnergy, p, mat, regi
204 } else {
205 res = GetRangeFromRestricteDEDX(kinEnergy,
206 }
207 return res;
208 }
209
210 //....oooOO0OOooo........oooOO0OOooo........oo
211
212 G4double G4EmCalculator::GetKinEnergy(G4double
213 const G4
214 const G4
215 const G4
216 {
217 G4double res = 0.0;
218 const G4MaterialCutsCouple* couple = FindCou
219 if(nullptr != couple && UpdateParticle(p, 1.
220 res = manager->GetEnergy(p, range, couple)
221 if(verbose>0) {
222 G4cout << "G4EmCalculator::GetKinEnergy:
223 << " KinE(MeV)= " << res/MeV
224 << " " << p->GetParticleName()
225 << " in " << mat->GetName()
226 << G4endl;
227 }
228 }
229 return res;
230 }
231
232 //....oooOO0OOooo........oooOO0OOooo........oo
233
234 G4double G4EmCalculator::GetCrossSectionPerVol
235 co
236 co
237 co
238 co
239 {
240 G4double res = 0.0;
241 const G4MaterialCutsCouple* couple = FindCou
242
243 if(nullptr != couple && UpdateParticle(p, ki
244 if(FindEmModel(p, processName, kinEnergy))
245 G4int idx = couple->GetIndex();
246 G4int procType = -1;
247 FindLambdaTable(p, processName, kinEnerg
248
249 G4VEmProcess* emproc = FindDiscreteProce
250 if(nullptr != emproc) {
251 res = emproc->GetCrossSection(kinEnergy, cou
252 } else if(currentLambda) {
253 // special tables are built for Msc mo
254 // procType is set in FindLambdaTable
255 if(procType==2) {
256 auto mscM = static_cast<G4VMscModel*
257 mscM->SetCurrentCouple(couple);
258 G4double tr1Mfp = mscM->GetTransport
259 if (tr1Mfp<DBL_MAX) {
260 res = 1./tr1Mfp;
261 }
262 } else {
263 G4double e = kinEnergy*massRatio;
264 res = (((*currentLambda)[idx])->Valu
265 }
266 } else {
267 res = ComputeCrossSectionPerVolume(kin
268 }
269 if(verbose>0) {
270 G4cout << "G4EmCalculator::GetXSPerVol
271 << " cross(cm-1)= " << res*cm
272 << " " << p->GetParticleName(
273 << " in " << mat->GetName();
274 if(verbose>1)
275 G4cout << " idx= " << idx << " Esc
276 << kinEnergy*massRatio
277 << " q2= " << chargeSquare;
278 G4cout << G4endl;
279 }
280 }
281 }
282 return res;
283 }
284
285 //....oooOO0OOooo........oooOO0OOooo........oo
286
287 G4double G4EmCalculator::GetShellIonisationCro
288 const
289 G4int
290 G4Ato
291 G4dou
292 {
293 G4double res = 0.0;
294 const G4ParticleDefinition* p = FindParticle
295 G4VAtomDeexcitation* ad = manager->AtomDeexc
296 if(nullptr != p && nullptr != ad) {
297 res = ad->GetShellIonisationCrossSectionPe
298 }
299 return res;
300 }
301
302 //....oooOO0OOooo........oooOO0OOooo........oo
303
304 G4double G4EmCalculator::GetMeanFreePath(G4dou
305 const
306 const
307 const
308 const
309 {
310 G4double res = DBL_MAX;
311 G4double x = GetCrossSectionPerVolume(kinEne
312 if(x > 0.0) { res = 1.0/x; }
313 if(verbose>1) {
314 G4cout << "G4EmCalculator::GetMeanFreePath
315 << " MFP(mm)= " << res/mm
316 << " " << p->GetParticleName()
317 << " in " << mat->GetName()
318 << G4endl;
319 }
320 return res;
321 }
322
323 //....oooOO0OOooo........oooOO0OOooo........oo
324
325 void G4EmCalculator::PrintDEDXTable(const G4Pa
326 {
327 const G4VEnergyLossProcess* elp = manager->G
328 G4cout << "##### DEDX Table for " << p->GetP
329 if(nullptr != elp) G4cout << *(elp->DEDXTabl
330 }
331
332 //....oooOO0OOooo........oooOO0OOooo........oo
333
334 void G4EmCalculator::PrintRangeTable(const G4P
335 {
336 const G4VEnergyLossProcess* elp = manager->G
337 G4cout << "##### Range Table for " << p->Get
338 if(nullptr != elp) G4cout << *(elp->RangeTab
339 }
340
341 //....oooOO0OOooo........oooOO0OOooo........oo
342
343 void G4EmCalculator::PrintInverseRangeTable(co
344 {
345 const G4VEnergyLossProcess* elp = manager->G
346 G4cout << "### G4EmCalculator: Inverse Range
347 << p->GetParticleName() << G4endl;
348 if(nullptr != elp) G4cout << *(elp->InverseR
349 }
350
351 //....oooOO0OOooo........oooOO0OOooo........oo
352
353 G4double G4EmCalculator::ComputeDEDX(G4double
354 const G4P
355 const G4S
356 const G4M
357 G4d
358 {
359 SetupMaterial(mat);
360 G4double res = 0.0;
361 if(verbose > 1) {
362 G4cout << "### G4EmCalculator::ComputeDEDX
363 << " in " << currentMaterialName
364 << " e(MeV)= " << kinEnergy/MeV <<
365 << G4endl;
366 }
367 if(UpdateParticle(p, kinEnergy)) {
368 if(FindEmModel(p, processName, kinEnergy))
369 G4double escaled = kinEnergy*massRatio;
370 if(nullptr != baseParticle) {
371 res = currentModel->ComputeDEDXPerVolume(mat
372
373 if(verbose > 1) {
374 G4cout << "Particle: " << p->GetParticleNa
375 << " E(MeV)=" << kinEnergy
376 << " Base particle: " << baseParticle->Ge
377 << " Escaled(MeV)= " << escaled
378 << " q2=" << chargeSquare << G4endl;
379 }
380 } else {
381 res = currentModel->ComputeDEDXPerVolume(mat
382 if(verbose > 1) {
383 G4cout << "Particle: " << p->GetParticleNa
384 << " E(MeV)=" << kinEnergy << G4endl;
385 }
386 }
387 if(verbose > 1) {
388 G4cout << currentModel->GetName() << ": DEDX
389 << " DEDX(MeV*cm^2/g)= "
390 << res*gram/(MeV*cm2*mat->GetDensity(
391 << G4endl;
392 }
393 // emulate smoothing procedure
394 if(applySmoothing && nullptr != loweMode
395 G4double eth = currentModel->LowEnergyLimit(
396 G4double res0 = 0.0;
397 G4double res1 = 0.0;
398 if(nullptr != baseParticle) {
399 res1 = chargeSquare*
400 currentModel->ComputeDEDXPerVolume(mat,
401 res0 = chargeSquare*
402 loweModel->ComputeDEDXPerVolume(mat, bas
403 } else {
404 res1 = currentModel->ComputeDEDXPerVolume(
405 res0 = loweModel->ComputeDEDXPerVolume(mat
406 }
407 if(res1 > 0.0 && escaled > 0.0) {
408 res *= (1.0 + (res0/res1 - 1.0)*eth/escale
409 }
410 if(verbose > 1) {
411 G4cout << "At boundary energy(MeV)= " << e
412 << " DEDX(MeV/mm)= " << res0*mm/MeV << "
413 << " after correction DEDX(MeV/mm)=" << r
414 }
415 }
416 // correction for ions
417 if(isIon) {
418 const G4double length = CLHEP::nm;
419 if(UpdateCouple(mat, cut)) {
420 G4double eloss = res*length;
421 dynParticle->SetKineticEnergy(kinEnergy);
422 currentModel->CorrectionsAlongStep(current
423 l
424 res = eloss/length;
425
426 if(verbose > 1) {
427 G4cout << "After Corrections: DEDX(MeV/m
428 << " DEDX(MeV*cm^2/g)= "
429 << res*gram/(MeV*cm2*mat->GetDensity())
430 }
431 }
432 }
433 if(verbose > 0) {
434 G4cout << "## E(MeV)= " << kinEnergy/MeV
435 << " DEDX(MeV/mm)= " << res*mm/MeV
436 << " DEDX(MeV*cm^2/g)= " << res*gram/
437 << " cut(MeV)= " << cut/MeV
438 << " " << p->GetParticleName()
439 << " in " << currentMaterialName
440 << " Zi^2= " << chargeSquare
441 << " isIon=" << isIon
442 << G4endl;
443 }
444 }
445 }
446 return res;
447 }
448
449 //....oooOO0OOooo........oooOO0OOooo........oo
450
451 G4double G4EmCalculator::ComputeElectronicDEDX
452
453
454
455 {
456 SetupMaterial(mat);
457 G4double dedx = 0.0;
458 if(UpdateParticle(part, kinEnergy)) {
459
460 G4LossTableManager* lManager = G4LossTable
461 const std::vector<G4VEnergyLossProcess*> v
462 lManager->GetEnergyLossProcessVector();
463 std::size_t n = vel.size();
464
465 //G4cout << "ComputeElectronicDEDX for " <
466 // << " n= " << n << G4endl;
467
468 for(std::size_t i=0; i<n; ++i) {
469 if(vel[i]) {
470 auto p = static_cast<G4VProcess*>(vel[
471 if(ActiveForParticle(part, p)) {
472 //G4cout << "idx= " << i << " " << (
473 // << " " << (vel[i])->Particle()-
474 dedx += ComputeDEDX(kinEnergy,part,(
475 }
476 }
477 }
478 }
479 return dedx;
480 }
481
482 //....oooOO0OOooo........oooOO0OOooo........oo
483
484 G4double
485 G4EmCalculator::ComputeDEDXForCutInRange(G4dou
486 const
487 const
488 G4dou
489 {
490 SetupMaterial(mat);
491 G4double dedx = 0.0;
492 if(UpdateParticle(part, kinEnergy)) {
493
494 G4LossTableManager* lManager = G4LossTable
495 const std::vector<G4VEnergyLossProcess*> v
496 lManager->GetEnergyLossProcessVector();
497 std::size_t n = vel.size();
498
499 if(mat != cutMaterial) {
500 cutMaterial = mat;
501 cutenergy[0] =
502 ComputeEnergyCutFromRangeCut(rangecut,
503 cutenergy[1] =
504 ComputeEnergyCutFromRangeCut(rangecut,
505 cutenergy[2] =
506 ComputeEnergyCutFromRangeCut(rangecut,
507 }
508
509 //G4cout << "ComputeElectronicDEDX for " <
510 // << " n= " << n << G4endl;
511
512 for(std::size_t i=0; i<n; ++i) {
513 if(vel[i]) {
514 auto p = static_cast<G4VProcess*>(vel[
515 if(ActiveForParticle(part, p)) {
516 //G4cout << "idx= " << i << " " << (
517 // << " " << (vel[i])->Particle()->
518 const G4ParticleDefinition* sec = (v
519 std::size_t idx = 0;
520 if(sec == G4Electron::Electron()) {
521 else if(sec == G4Positron::Positron(
522
523 dedx += ComputeDEDX(kinEnergy,part,(
524 mat,cutenergy[id
525 }
526 }
527 }
528 }
529 return dedx;
530 }
531
532 //....oooOO0OOooo........oooOO0OOooo........oo
533
534 G4double G4EmCalculator::ComputeTotalDEDX(G4do
535 cons
536 cons
537 G4do
538 {
539 G4double dedx = ComputeElectronicDEDX(kinEne
540 if(mass > 700.*MeV) { dedx += ComputeNuclear
541 return dedx;
542 }
543
544 //....oooOO0OOooo........oooOO0OOooo........oo
545
546 G4double G4EmCalculator::ComputeNuclearDEDX(G4
547 const G4
548 const G4
549 {
550 G4double res = 0.0;
551 G4VEmProcess* nucst = FindDiscreteProcess(p,
552 if(nucst) {
553 G4VEmModel* mod = nucst->EmModel();
554 if(mod) {
555 mod->SetFluctuationFlag(false);
556 res = mod->ComputeDEDXPerVolume(mat, p,
557 }
558 }
559
560 if(verbose > 1) {
561 G4cout << p->GetParticleName() << " E(MeV
562 << " NuclearDEDX(MeV/mm)= " << res*
563 << " NuclearDEDX(MeV*cm^2/g)= "
564 << res*gram/(MeV*cm2*mat->GetDensit
565 << G4endl;
566 }
567 return res;
568 }
569
570 //....oooOO0OOooo........oooOO0OOooo........oo
571
572 G4double G4EmCalculator::ComputeCrossSectionPe
573
574 c
575 c
576 c
577
578 {
579 SetupMaterial(mat);
580 G4double res = 0.0;
581 if(UpdateParticle(p, kinEnergy)) {
582 if(FindEmModel(p, processName, kinEnergy))
583 G4double e = kinEnergy;
584 G4double aCut = std::max(cut, theParamet
585 if(baseParticle) {
586 e *= kinEnergy*massRatio;
587 res = currentModel->CrossSectionPerVol
588 mat, baseParticle, e, aCut, e) *
589 } else {
590 res = currentModel->CrossSectionPerVol
591 }
592 if(verbose>0) {
593 G4cout << "G4EmCalculator::ComputeXSPe
594 << kinEnergy/MeV
595 << " cross(cm-1)= " << res*cm
596 << " cut(keV)= " << aCut/keV
597 << " " << p->GetParticleName(
598 << " in " << mat->GetName()
599 << G4endl;
600 }
601 }
602 }
603 return res;
604 }
605
606 //....oooOO0OOooo........oooOO0OOooo........oo
607
608 G4double
609 G4EmCalculator::ComputeCrossSectionPerAtom(G4d
610 con
611 con
612 G4d
613 G4d
614 {
615 G4double res = 0.0;
616 if(UpdateParticle(p, kinEnergy)) {
617 G4int iz = G4lrint(Z);
618 CheckMaterial(iz);
619 if(FindEmModel(p, processName, kinEnergy))
620 G4double e = kinEnergy;
621 G4double aCut = std::max(cut, theParamet
622 if(baseParticle) {
623 e *= kinEnergy*massRatio;
624 currentModel->InitialiseForElement(bas
625 res = currentModel->ComputeCrossSectio
626 baseParticle, e, Z, A, aCut) * c
627 } else {
628 currentModel->InitialiseForElement(p,
629 res = currentModel->ComputeCrossSectio
630 }
631 if(verbose>0) {
632 G4cout << "E(MeV)= " << kinEnergy/MeV
633 << " cross(barn)= " << res/barn
634 << " " << p->GetParticleName(
635 << " Z= " << Z << " A= " << A/
636 << " cut(keV)= " << aCut/keV
637 << G4endl;
638 }
639 }
640 }
641 return res;
642 }
643
644 //....oooOO0OOooo........oooOO0OOooo........oo
645
646 G4double
647 G4EmCalculator::ComputeCrossSectionPerShell(G4
648 co
649 co
650 G4
651 G4
652 {
653 G4double res = 0.0;
654 if(UpdateParticle(p, kinEnergy)) {
655 CheckMaterial(Z);
656 if(FindEmModel(p, processName, kinEnergy))
657 G4double e = kinEnergy;
658 G4double aCut = std::max(cut, theParamet
659 if(nullptr != baseParticle) {
660 e *= kinEnergy*massRatio;
661 currentModel->InitialiseForElement(bas
662 res =
663 currentModel->ComputeCrossSectionPer
664
665 } else {
666 currentModel->InitialiseForElement(p,
667 res = currentModel->ComputeCrossSectio
668 }
669 if(verbose>0) {
670 G4cout << "E(MeV)= " << kinEnergy/MeV
671 << " cross(barn)= " << res/barn
672 << " " << p->GetParticleName(
673 << " Z= " << Z << " shellIdx=
674 << " cut(keV)= " << aCut/keV
675 << G4endl;
676 }
677 }
678 }
679 return res;
680 }
681
682 //....oooOO0OOooo........oooOO0OOooo........oo
683
684 G4double
685 G4EmCalculator::ComputeGammaAttenuationLength(
686
687 {
688 G4double res = 0.0;
689 const G4ParticleDefinition* gamma = G4Gamma:
690 res += ComputeCrossSectionPerVolume(kinEnerg
691 res += ComputeCrossSectionPerVolume(kinEnerg
692 res += ComputeCrossSectionPerVolume(kinEnerg
693 res += ComputeCrossSectionPerVolume(kinEnerg
694 if(res > 0.0) { res = 1.0/res; }
695 return res;
696 }
697
698 //....oooOO0OOooo........oooOO0OOooo........oo
699
700 G4double G4EmCalculator::ComputeShellIonisatio
701 const
702 G4int
703 G4Ato
704 G4dou
705 const
706 {
707 G4double res = 0.0;
708 const G4ParticleDefinition* p = FindParticle
709 G4VAtomDeexcitation* ad = manager->AtomDeexc
710 if(p && ad) {
711 res = ad->ComputeShellIonisationCrossSecti
712
713 }
714 return res;
715 }
716
717 //....oooOO0OOooo........oooOO0OOooo........oo
718
719 G4double G4EmCalculator::ComputeMeanFreePath(G
720 c
721 c
722 c
723 G
724 {
725 G4double mfp = DBL_MAX;
726 G4double x =
727 ComputeCrossSectionPerVolume(kinEnergy, p,
728 if(x > 0.0) { mfp = 1.0/x; }
729 if(verbose>1) {
730 G4cout << "E(MeV)= " << kinEnergy/MeV
731 << " MFP(mm)= " << mfp/mm
732 << " " << p->GetParticleName()
733 << " in " << mat->GetName()
734 << G4endl;
735 }
736 return mfp;
737 }
738
739 //....oooOO0OOooo........oooOO0OOooo........oo
740
741 G4double G4EmCalculator::ComputeEnergyCutFromR
742 G4double range,
743 const G4ParticleDefin
744 const G4Material* mat
745 {
746 return G4ProductionCutsTable::GetProductionC
747 ConvertRangeToEnergy(part, mat, range);
748 }
749
750 //....oooOO0OOooo........oooOO0OOooo........oo
751
752 G4bool G4EmCalculator::UpdateParticle(const G4
753 G4double
754 {
755 if(p != currentParticle) {
756
757 // new particle
758 currentParticle = p;
759 dynParticle->SetDefinition(const_cast<G4Pa
760 dynParticle->SetKineticEnergy(kinEnergy);
761 baseParticle = nullptr;
762 currentParticleName = p->GetParticleName()
763 massRatio = 1.0;
764 mass = p->GetPDGMass();
765 chargeSquare = 1.0;
766 currentProcess = manager->GetEnergyLossPr
767 currentProcessName = "";
768 isIon = false;
769
770 // ionisation process exist
771 if(nullptr != currentProcess) {
772 currentProcessName = currentProcess->Get
773 baseParticle = currentProcess->BaseParti
774 if(currentProcessName == "ionIoni" && p-
775 baseParticle = theGenericIon;
776 isIon = true;
777 }
778
779 // base particle is used
780 if(nullptr != baseParticle) {
781 massRatio = baseParticle->GetPDGMass()
782 G4double q = p->GetPDGCharge()/basePar
783 chargeSquare = q*q;
784 }
785 }
786 }
787 // Effective charge for ions
788 if(isIon && nullptr != currentProcess) {
789 chargeSquare =
790 corr->EffectiveChargeSquareRatio(p, curr
791 currentProcess->SetDynamicMassCharge(massR
792 if(verbose>1) {
793 G4cout <<"\n NewIon: massR= "<< massRati
794 << chargeSquare << " " << currentProce
795 }
796 }
797 return true;
798 }
799
800 //....oooOO0OOooo........oooOO0OOooo........oo
801
802 const G4ParticleDefinition* G4EmCalculator::Fi
803 {
804 const G4ParticleDefinition* p = nullptr;
805 if(name != currentParticleName) {
806 p = G4ParticleTable::GetParticleTable()->F
807 if(nullptr == p) {
808 G4cout << "### WARNING: G4EmCalculator::
809 << name << G4endl;
810 }
811 } else {
812 p = currentParticle;
813 }
814 return p;
815 }
816
817 //....oooOO0OOooo........oooOO0OOooo........oo
818
819 const G4ParticleDefinition* G4EmCalculator::Fi
820 {
821 const G4ParticleDefinition* p = ionTable->Ge
822 return p;
823 }
824
825 //....oooOO0OOooo........oooOO0OOooo........oo
826
827 const G4Material* G4EmCalculator::FindMaterial
828 {
829 if(name != currentMaterialName) {
830 SetupMaterial(G4Material::GetMaterial(name
831 if(nullptr == currentMaterial) {
832 G4cout << "### WARNING: G4EmCalculator::
833 << name << G4endl;
834 }
835 }
836 return currentMaterial;
837 }
838
839 //....oooOO0OOooo........oooOO0OOooo........oo
840
841 const G4Region* G4EmCalculator::FindRegion(con
842 {
843 return G4EmUtility::FindRegion(reg);
844 }
845
846 //....oooOO0OOooo........oooOO0OOooo........oo
847
848 const G4MaterialCutsCouple* G4EmCalculator::Fi
849 const G4Material*
850 const G4Region* re
851 {
852 const G4MaterialCutsCouple* couple = nullptr
853 SetupMaterial(material);
854 if(nullptr != currentMaterial) {
855 // Access to materials
856 const G4ProductionCutsTable* theCoupleTabl
857 G4ProductionCutsTable::GetProductionCuts
858 const G4Region* r = region;
859 if(nullptr != r) {
860 couple = theCoupleTable->GetMaterialCuts
861
862 } else {
863 G4RegionStore* store = G4RegionStore::Ge
864 std::size_t nr = store->size();
865 if(0 < nr) {
866 for(std::size_t i=0; i<nr; ++i) {
867 couple = theCoupleTable->GetMaterial
868 material, ((*store)[i])->GetProduc
869 if(nullptr != couple) { break; }
870 }
871 }
872 }
873 }
874 if(nullptr == couple) {
875 G4ExceptionDescription ed;
876 ed << "G4EmCalculator::FindCouple: fail fo
877 << currentMaterialName << ">";
878 if(region) { ed << " and region " << regio
879 G4Exception("G4EmCalculator::FindCouple",
880 FatalException, ed);
881 }
882 return couple;
883 }
884
885 //....oooOO0OOooo........oooOO0OOooo........oo
886
887 G4bool G4EmCalculator::UpdateCouple(const G4Ma
888 {
889 SetupMaterial(material);
890 if(!currentMaterial) { return false; }
891 for (G4int i=0; i<nLocalMaterials; ++i) {
892 if(material == localMaterials[i] && cut ==
893 currentCouple = localCouples[i];
894 currentCoupleIndex = currentCouple->GetI
895 currentCut = cut;
896 return true;
897 }
898 }
899 const G4MaterialCutsCouple* cc = new G4Mater
900 localMaterials.push_back(material);
901 localCouples.push_back(cc);
902 localCuts.push_back(cut);
903 ++nLocalMaterials;
904 currentCouple = cc;
905 currentCoupleIndex = currentCouple->GetIndex
906 currentCut = cut;
907 return true;
908 }
909
910 //....oooOO0OOooo........oooOO0OOooo........oo
911
912 void G4EmCalculator::FindLambdaTable(const G4P
913 const G4S
914 G4double
915 {
916 // Search for the process
917 if (!currentLambda || p != lambdaParticle ||
918 lambdaName = processName;
919 currentLambda = nullptr;
920 lambdaParticle = p;
921 isApplicable = false;
922
923 const G4ParticleDefinition* part = (isIon)
924
925 // Search for energy loss process
926 currentName = processName;
927 currentModel = nullptr;
928 loweModel = nullptr;
929
930 G4VEnergyLossProcess* elproc = FindEnLossP
931 if(nullptr != elproc) {
932 currentLambda = elproc->LambdaTable();
933 proctype = 0;
934 if(nullptr != currentLambda) {
935 isApplicable = true;
936 if(verbose>1) {
937 G4cout << "G4VEnergyLossProcess is f
938 << G4endl;
939 }
940 }
941 curProcess = elproc;
942 return;
943 }
944
945 // Search for discrete process
946 G4VEmProcess* proc = FindDiscreteProcess(p
947 if(nullptr != proc) {
948 currentLambda = proc->LambdaTable();
949 proctype = 1;
950 if(nullptr != currentLambda) {
951 isApplicable = true;
952 if(verbose>1) {
953 G4cout << "G4VEmProcess is found out
954 }
955 }
956 curProcess = proc;
957 return;
958 }
959
960 // Search for msc process
961 G4VMultipleScattering* msc = FindMscProces
962 if(nullptr != msc) {
963 currentModel = msc->SelectModel(kinEnerg
964 proctype = 2;
965 if(nullptr != currentModel) {
966 currentLambda = currentModel->GetCross
967 if(nullptr != currentLambda) {
968 isApplicable = true;
969 if(verbose>1) {
970 G4cout << "G4VMultipleScattering i
971 << G4endl;
972 }
973 }
974 }
975 curProcess = msc;
976 }
977 }
978 }
979
980 //....oooOO0OOooo........oooOO0OOooo........oo
981
982 G4bool G4EmCalculator::FindEmModel(const G4Par
983 const G4Str
984 G4
985 {
986 isApplicable = false;
987 if(nullptr == p || nullptr == currentMateria
988 G4cout << "G4EmCalculator::FindEmModel WAR
989 << " or materail defined; particle:
990 return isApplicable;
991 }
992 G4String partname = p->GetParticleName();
993 G4double scaledEnergy = kinEnergy*massRatio;
994 const G4ParticleDefinition* part = (isIon) ?
995
996 if(verbose > 1) {
997 G4cout << "## G4EmCalculator::FindEmModel
998 << " (type= " << p->GetParticleType
999 << ") and " << processName << " at
1000 << G4endl;
1001 if(p != part) { G4cout << " GenericIon i
1002 }
1003
1004 // Search for energy loss process
1005 currentName = processName;
1006 currentModel = nullptr;
1007 loweModel = nullptr;
1008 std::size_t idx = 0;
1009
1010 G4VEnergyLossProcess* elproc = FindEnLossPr
1011 if(nullptr != elproc) {
1012 currentModel = elproc->SelectModelForMate
1013 currentModel->InitialiseForMaterial(part,
1014 currentModel->SetupForMaterial(part, curr
1015 G4double eth = currentModel->LowEnergyLim
1016 if(eth > 0.0) {
1017 loweModel = elproc->SelectModelForMater
1018 if(loweModel == currentModel) { loweMod
1019 else {
1020 loweModel->InitialiseForMaterial(part
1021 loweModel->SetupForMaterial(part, cur
1022 }
1023 }
1024 }
1025
1026 // Search for discrete process
1027 if(nullptr == currentModel) {
1028 G4VEmProcess* proc = FindDiscreteProcess(
1029 if(nullptr != proc) {
1030 currentModel = proc->SelectModelForMate
1031 currentModel->InitialiseForMaterial(par
1032 currentModel->SetupForMaterial(part, cu
1033 G4double eth = currentModel->LowEnergyL
1034 if(eth > 0.0) {
1035 loweModel = proc->SelectModelForMater
1036 if(loweModel == currentModel) { loweM
1037 else {
1038 loweModel->InitialiseForMaterial(pa
1039 loweModel->SetupForMaterial(part, c
1040 }
1041 }
1042 }
1043 }
1044
1045 // Search for msc process
1046 if(nullptr == currentModel) {
1047 G4VMultipleScattering* proc = FindMscProc
1048 if(nullptr != proc) {
1049 currentModel = proc->SelectModel(kinEne
1050 loweModel = nullptr;
1051 }
1052 }
1053 if(nullptr != currentModel) {
1054 if(loweModel == currentModel) { loweModel
1055 isApplicable = true;
1056 currentModel->InitialiseForMaterial(part,
1057 if(loweModel) {
1058 loweModel->InitialiseForMaterial(part,
1059 }
1060 if(verbose > 1) {
1061 G4cout << " Model <" << currentModel-
1062 << "> Emin(MeV)= " << currentMod
1063 << " for " << part->GetParticleN
1064 if(nullptr != elproc) {
1065 G4cout << " and " << elproc->GetProce
1066 << G4endl;
1067 }
1068 if(nullptr != loweModel) {
1069 G4cout << " LowEnergy model <" << low
1070 }
1071 G4cout << G4endl;
1072 }
1073 }
1074 return isApplicable;
1075 }
1076
1077 //....oooOO0OOooo........oooOO0OOooo........o
1078
1079 G4VEnergyLossProcess*
1080 G4EmCalculator::FindEnLossProcess(const G4Par
1081 const G4Str
1082 {
1083 G4VEnergyLossProcess* proc = nullptr;
1084 const std::vector<G4VEnergyLossProcess*> v
1085 manager->GetEnergyLossProcessVector();
1086 std::size_t n = v.size();
1087 for(std::size_t i=0; i<n; ++i) {
1088 if((v[i])->GetProcessName() == processNam
1089 auto p = static_cast<G4VProcess*>(v[i])
1090 if(ActiveForParticle(part, p)) {
1091 proc = v[i];
1092 break;
1093 }
1094 }
1095 }
1096 return proc;
1097 }
1098
1099 //....oooOO0OOooo........oooOO0OOooo........o
1100
1101 G4VEmProcess*
1102 G4EmCalculator::FindDiscreteProcess(const G4P
1103 const G4S
1104 {
1105 G4VEmProcess* proc = nullptr;
1106 auto v = manager->GetEmProcessVector();
1107 std::size_t n = v.size();
1108 for(std::size_t i=0; i<n; ++i) {
1109 const G4String& pName = v[i]->GetProcessN
1110 if(pName == "GammaGeneralProc") {
1111 proc = v[i]->GetEmProcess(processName);
1112 break;
1113 } else if(pName == processName) {
1114 const auto p = static_cast<G4VProcess*>
1115 if(ActiveForParticle(part, p)) {
1116 proc = v[i];
1117 break;
1118 }
1119 }
1120 }
1121 return proc;
1122 }
1123
1124 //....oooOO0OOooo........oooOO0OOooo........o
1125
1126 G4VMultipleScattering*
1127 G4EmCalculator::FindMscProcess(const G4Partic
1128 const G4String
1129 {
1130 G4VMultipleScattering* proc = nullptr;
1131 const std::vector<G4VMultipleScattering*> v
1132 manager->GetMultipleScatteringVector();
1133 std::size_t n = v.size();
1134 for(std::size_t i=0; i<n; ++i) {
1135 if((v[i])->GetProcessName() == processNam
1136 auto p = static_cast<G4VProcess*>(v[i])
1137 if(ActiveForParticle(part, p)) {
1138 proc = v[i];
1139 break;
1140 }
1141 }
1142 }
1143 return proc;
1144 }
1145
1146 //....oooOO0OOooo........oooOO0OOooo........o
1147
1148 G4VProcess* G4EmCalculator::FindProcess(const
1149 const
1150 {
1151 G4VProcess* proc = nullptr;
1152 const G4ProcessManager* procman = part->Get
1153 G4ProcessVector* pv = procman->GetProcessLi
1154 G4int nproc = (G4int)pv->size();
1155 for(G4int i=0; i<nproc; ++i) {
1156 if(processName == (*pv)[i]->GetProcessNam
1157 proc = (*pv)[i];
1158 break;
1159 }
1160 }
1161 return proc;
1162 }
1163
1164 //....oooOO0OOooo........oooOO0OOooo........o
1165
1166 G4bool G4EmCalculator::ActiveForParticle(cons
1167 G4VP
1168 {
1169 G4ProcessManager* pm = part->GetProcessMana
1170 G4ProcessVector* pv = pm->GetProcessList();
1171 G4int n = (G4int)pv->size();
1172 G4bool res = false;
1173 for(G4int i=0; i<n; ++i) {
1174 if((*pv)[i] == proc) {
1175 if(pm->GetProcessActivation(i)) { res =
1176 break;
1177 }
1178 }
1179 return res;
1180 }
1181
1182 //....oooOO0OOooo........oooOO0OOooo........o
1183
1184 void G4EmCalculator::SetupMaterial(const G4Ma
1185 {
1186 if(mat) {
1187 currentMaterial = mat;
1188 currentMaterialName = mat->GetName();
1189 } else {
1190 currentMaterial = nullptr;
1191 currentMaterialName = "";
1192 }
1193 }
1194
1195 //....oooOO0OOooo........oooOO0OOooo........o
1196
1197 void G4EmCalculator::SetupMaterial(const G4St
1198 {
1199 SetupMaterial(nist->FindOrBuildMaterial(mna
1200 }
1201
1202 //....oooOO0OOooo........oooOO0OOooo........o
1203
1204 void G4EmCalculator::CheckMaterial(G4int Z)
1205 {
1206 G4bool isFound = false;
1207 if(nullptr != currentMaterial) {
1208 G4int nn = (G4int)currentMaterial->GetNum
1209 for(G4int i=0; i<nn; ++i) {
1210 if(Z == currentMaterial->GetElement(i)-
1211 isFound = true;
1212 break;
1213 }
1214 }
1215 }
1216 if(!isFound) {
1217 SetupMaterial(nist->FindOrBuildSimpleMate
1218 }
1219 }
1220
1221 //....oooOO0OOooo........oooOO0OOooo........o
1222
1223 void G4EmCalculator::SetVerbose(G4int verb)
1224 {
1225 verbose = verb;
1226 }
1227
1228 //....oooOO0OOooo........oooOO0OOooo........o
1229
1230