Geant4 10.7.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4EmCalculator.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// -------------------------------------------------------------------
28//
29// GEANT4 Class 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........oooOO0OOooo........oooOO0OOooo....
44//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46#include "G4EmCalculator.hh"
47#include "G4SystemOfUnits.hh"
48#include "G4LossTableManager.hh"
49#include "G4EmParameters.hh"
50#include "G4NistManager.hh"
51#include "G4VEmProcess.hh"
52#include "G4VEnergyLossProcess.hh"
53#include "G4VMultipleScattering.hh"
54#include "G4Material.hh"
55#include "G4MaterialCutsCouple.hh"
56#include "G4ParticleDefinition.hh"
57#include "G4ParticleTable.hh"
58#include "G4IonTable.hh"
59#include "G4PhysicsTable.hh"
60#include "G4ProductionCutsTable.hh"
61#include "G4ProcessManager.hh"
62#include "G4ionEffectiveCharge.hh"
63#include "G4RegionStore.hh"
64#include "G4Element.hh"
65#include "G4EmCorrections.hh"
66#include "G4GenericIon.hh"
67#include "G4ProcessVector.hh"
68#include "G4Gamma.hh"
69#include "G4Electron.hh"
70#include "G4Positron.hh"
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
74G4EmCalculator::G4EmCalculator()
75{
76 manager = G4LossTableManager::Instance();
77 nist = G4NistManager::Instance();
78 theParameters = G4EmParameters::Instance();
79 corr = manager->EmCorrections();
80 nLocalMaterials = 0;
81 verbose = 0;
82 currentCoupleIndex = 0;
83 currentCouple = nullptr;
84 currentMaterial = cutMaterial = nullptr;
85 currentParticle = nullptr;
86 lambdaParticle = nullptr;
87 baseParticle = nullptr;
88 currentLambda = nullptr;
89 currentModel = nullptr;
90 currentProcess = nullptr;
91 curProcess = nullptr;
92 loweModel = nullptr;
93 chargeSquare = 1.0;
94 massRatio = 1.0;
95 mass = 0.0;
96 currentCut = 0.0;
97 cutenergy[0] = cutenergy[1] = cutenergy[2] = DBL_MAX;
98 currentParticleName= "";
99 currentMaterialName= "";
100 currentName = "";
101 lambdaName = "";
102 theGenericIon = G4GenericIon::GenericIon();
103 ionEffCharge = new G4ionEffectiveCharge();
104 ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
105 isIon = false;
106 isApplicable = false;
107}
108
109//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
110
111G4EmCalculator::~G4EmCalculator()
112{
113 delete ionEffCharge;
114 for (G4int i=0; i<nLocalMaterials; ++i) {
115 delete localCouples[i];
116 }
117}
118
119//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
120
121G4double G4EmCalculator::GetDEDX(G4double kinEnergy,
122 const G4ParticleDefinition* p,
123 const G4Material* mat,
124 const G4Region* region)
125{
126 G4double res = 0.0;
127 const G4MaterialCutsCouple* couple = FindCouple(mat, region);
128 if(couple && UpdateParticle(p, kinEnergy) ) {
129 res = manager->GetDEDX(p, kinEnergy, couple);
130
131 if(isIon) {
132 if(FindEmModel(p, currentProcessName, kinEnergy)) {
133 G4double length = CLHEP::nm;
134 G4double eloss = res*length;
135 //G4cout << "### GetDEDX: E= " << kinEnergy << " dedx0= " << res
136 // << " de= " << eloss << G4endl;;
137 G4double niel = 0.0;
138 dynParticle.SetKineticEnergy(kinEnergy);
139 currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
140 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
141 res = eloss/length;
142 //G4cout << " de1= " << eloss << " res1= " << res
143 // << " " << p->GetParticleName() <<G4endl;;
144 }
145 }
146
147 if(verbose>0) {
148 G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
149 << " DEDX(MeV/mm)= " << res*mm/MeV
150 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
151 << " " << p->GetParticleName()
152 << " in " << mat->GetName()
153 << " isIon= " << isIon
154 << G4endl;
155 }
156 }
157 return res;
158}
159
160//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161
162G4double G4EmCalculator::GetRangeFromRestricteDEDX(G4double kinEnergy,
163 const G4ParticleDefinition* p,
164 const G4Material* mat,
165 const G4Region* region)
166{
167 G4double res = 0.0;
168 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
169 if(couple && UpdateParticle(p, kinEnergy)) {
170 res = manager->GetRangeFromRestricteDEDX(p, kinEnergy, couple);
171 if(verbose>1) {
172 G4cout << " G4EmCalculator::GetRangeFromRestrictedDEDX: E(MeV)= "
173 << kinEnergy/MeV
174 << " range(mm)= " << res/mm
175 << " " << p->GetParticleName()
176 << " in " << mat->GetName()
177 << G4endl;
178 }
179 }
180 return res;
181}
182
183//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
184
185G4double G4EmCalculator::GetCSDARange(G4double kinEnergy,
186 const G4ParticleDefinition* p,
187 const G4Material* mat,
188 const G4Region* region)
189{
190 G4double res = 0.0;
191 if(!theParameters->BuildCSDARange()) {
192 G4ExceptionDescription ed;
193 ed << "G4EmCalculator::GetCSDARange: CSDA table is not built; "
194 << " use UI command: /process/eLoss/CSDARange true";
195 G4Exception("G4EmCalculator::GetCSDARange", "em0077",
196 JustWarning, ed);
197 return res;
198 }
199
200 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
201 if(couple && UpdateParticle(p, kinEnergy)) {
202 res = manager->GetCSDARange(p, kinEnergy, couple);
203 if(verbose>1) {
204 G4cout << " G4EmCalculator::GetCSDARange: E(MeV)= " << kinEnergy/MeV
205 << " range(mm)= " << res/mm
206 << " " << p->GetParticleName()
207 << " in " << mat->GetName()
208 << G4endl;
209 }
210 }
211 return res;
212}
213
214//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
215
216G4double G4EmCalculator::GetRange(G4double kinEnergy,
217 const G4ParticleDefinition* p,
218 const G4Material* mat,
219 const G4Region* region)
220{
221 G4double res = 0.0;
222 if(theParameters->BuildCSDARange()) {
223 res = GetCSDARange(kinEnergy, p, mat, region);
224 } else {
225 res = GetRangeFromRestricteDEDX(kinEnergy, p, mat, region);
226 }
227 return res;
228}
229
230//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
231
232G4double G4EmCalculator::GetKinEnergy(G4double range,
233 const G4ParticleDefinition* p,
234 const G4Material* mat,
235 const G4Region* region)
236{
237 G4double res = 0.0;
238 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
239 if(couple && UpdateParticle(p, 1.0*GeV)) {
240 res = manager->GetEnergy(p, range, couple);
241 if(verbose>0) {
242 G4cout << "G4EmCalculator::GetKinEnergy: Range(mm)= " << range/mm
243 << " KinE(MeV)= " << res/MeV
244 << " " << p->GetParticleName()
245 << " in " << mat->GetName()
246 << G4endl;
247 }
248 }
249 return res;
250}
251
252//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
253
254G4double G4EmCalculator::GetCrossSectionPerVolume(G4double kinEnergy,
255 const G4ParticleDefinition* p,
256 const G4String& processName,
257 const G4Material* mat,
258 const G4Region* region)
259{
260 G4double res = 0.0;
261 const G4MaterialCutsCouple* couple = FindCouple(mat,region);
262
263 if(couple && UpdateParticle(p, kinEnergy)) {
264 if(FindEmModel(p, processName, kinEnergy)) {
265 G4int idx = couple->GetIndex();
266 G4int procType = -1;
267 FindLambdaTable(p, processName, kinEnergy, procType);
268
269 G4VEmProcess* emproc = FindDiscreteProcess(p, processName);
270 if(emproc) {
271 res = emproc->CrossSectionPerVolume(kinEnergy, couple);
272 } else if(currentLambda) {
273 // special tables are built for Msc models (procType is set in FindLambdaTable
274 if(procType==2) {
275 G4VMscModel* mscM = static_cast<G4VMscModel*>(currentModel);
276 mscM->SetCurrentCouple(couple);
277 G4double tr1Mfp = mscM->GetTransportMeanFreePath(p, kinEnergy);
278 if (tr1Mfp<DBL_MAX) {
279 res = 1./tr1Mfp;
280 }
281 } else {
282 G4double e = kinEnergy*massRatio;
283 res = (((*currentLambda)[idx])->Value(e))*chargeSquare;
284 }
285 } else {
286 res = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, kinEnergy);
287 }
288 if(verbose>0) {
289 G4cout << "G4EmCalculator::GetXSPerVolume: E(MeV)= " << kinEnergy/MeV
290 << " cross(cm-1)= " << res*cm
291 << " " << p->GetParticleName()
292 << " in " << mat->GetName();
293 if(verbose>1)
294 G4cout << " idx= " << idx << " Escaled((MeV)= "
295 << kinEnergy*massRatio
296 << " q2= " << chargeSquare;
297 G4cout << G4endl;
298 }
299 }
300 }
301 return res;
302}
303
304//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
305
306G4double G4EmCalculator::GetShellIonisationCrossSectionPerAtom(
307 const G4String& particle,
308 G4int Z,
309 G4AtomicShellEnumerator shell,
310 G4double kinEnergy)
311{
312 G4double res = 0.0;
313 const G4ParticleDefinition* p = FindParticle(particle);
314 G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
315 if(p && ad) {
316 res = ad->GetShellIonisationCrossSectionPerAtom(p, Z, shell, kinEnergy);
317 }
318 return res;
319}
320
321//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
322
323G4double G4EmCalculator::GetMeanFreePath(G4double kinEnergy,
324 const G4ParticleDefinition* p,
325 const G4String& processName,
326 const G4Material* mat,
327 const G4Region* region)
328{
329 G4double res = DBL_MAX;
330 G4double x = GetCrossSectionPerVolume(kinEnergy,p, processName, mat,region);
331 if(x > 0.0) { res = 1.0/x; }
332 if(verbose>1) {
333 G4cout << "G4EmCalculator::GetMeanFreePath: E(MeV)= " << kinEnergy/MeV
334 << " MFP(mm)= " << res/mm
335 << " " << p->GetParticleName()
336 << " in " << mat->GetName()
337 << G4endl;
338 }
339 return res;
340}
341
342//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
343
344void G4EmCalculator::PrintDEDXTable(const G4ParticleDefinition* p)
345{
346 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
347 G4cout << "##### DEDX Table for " << p->GetParticleName() << G4endl;
348 if(elp) G4cout << *(elp->DEDXTable()) << G4endl;
349}
350
351//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352
353void G4EmCalculator::PrintRangeTable(const G4ParticleDefinition* p)
354{
355 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
356 G4cout << "##### Range Table for " << p->GetParticleName() << G4endl;
357 if(elp) G4cout << *(elp->RangeTableForLoss()) << G4endl;
358}
359
360//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
361
362void G4EmCalculator::PrintInverseRangeTable(const G4ParticleDefinition* p)
363{
364 const G4VEnergyLossProcess* elp = FindEnergyLossProcess(p);
365 G4cout << "### G4EmCalculator: Inverse Range Table for "
366 << p->GetParticleName() << G4endl;
367 if(elp) G4cout << *(elp->InverseRangeTable()) << G4endl;
368}
369
370//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
371
372G4double G4EmCalculator::ComputeDEDX(G4double kinEnergy,
373 const G4ParticleDefinition* p,
374 const G4String& processName,
375 const G4Material* mat,
376 G4double cut)
377{
378 SetupMaterial(mat);
379 G4double res = 0.0;
380 if(verbose > 1) {
381 G4cout << "### G4EmCalculator::ComputeDEDX: " << p->GetParticleName()
382 << " in " << currentMaterialName
383 << " e(MeV)= " << kinEnergy/MeV << " cut(MeV)= " << cut/MeV
384 << G4endl;
385 }
386 if(UpdateParticle(p, kinEnergy)) {
387 if(FindEmModel(p, processName, kinEnergy)) {
388
389 // Special case of ICRU'73 model
390 const G4String& mname = currentModel->GetName();
391 if(mname == "ParamICRU73" || mname == "LinhardSorensen" || mname == "Atima") {
392 res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
393 if(verbose > 1) {
394 G4cout << mname << " ion E(MeV)= " << kinEnergy << " ";
395 G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
396 << " DEDX(MeV*cm^2/g)= "
397 << res*gram/(MeV*cm2*mat->GetDensity())
398 << G4endl;
399 }
400 } else {
401
402 G4double escaled = kinEnergy*massRatio;
403 if(baseParticle) {
404 res = currentModel->ComputeDEDXPerVolume(
405 mat, baseParticle, escaled, cut) * chargeSquare;
406 if(verbose > 1) {
407 G4cout << baseParticle->GetParticleName()
408 << " Escaled(MeV)= " << escaled;
409 }
410 } else {
411 res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
412 if(verbose > 1) { G4cout << " no basePart E(MeV)= " << kinEnergy << " "; }
413 }
414 if(verbose > 1) {
415 G4cout << currentModel->GetName() << ": DEDX(MeV/mm)= " << res*mm/MeV
416 << " DEDX(MeV*cm^2/g)= "
417 << res*gram/(MeV*cm2*mat->GetDensity())
418 << G4endl;
419 }
420
421 // emulate smoothing procedure
422 G4double eth = currentModel->LowEnergyLimit();
423 // G4cout << "massRatio= " << massRatio << " eth= " << eth << G4endl;
424 if(loweModel) {
425 G4double res0 = 0.0;
426 G4double res1 = 0.0;
427 if(baseParticle) {
428 res1 = currentModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
429 * chargeSquare;
430 res0 = loweModel->ComputeDEDXPerVolume(mat, baseParticle, eth, cut)
431 * chargeSquare;
432 } else {
433 res1 = currentModel->ComputeDEDXPerVolume(mat, p, eth, cut);
434 res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
435 }
436 if(verbose > 1) {
437 G4cout << "At boundary energy(MeV)= " << eth/MeV
438 << " DEDX(MeV/mm)= " << res1*mm/MeV
439 << G4endl;
440 }
441
442 //G4cout << "eth= " << eth << " escaled= " << escaled
443 // << " res0= " << res0 << " res1= "
444 // << res1 << " q2= " << chargeSquare << G4endl;
445
446 if(res1 > 0.0 && escaled > 0.0) {
447 res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
448 }
449 }
450
451 // low energy correction for ions
452 if(isIon) {
453 G4double length = CLHEP::nm;
454 const G4Region* r = 0;
455 const G4MaterialCutsCouple* couple = FindCouple(mat, r);
456 G4double eloss = res*length;
457 G4double niel = 0.0;
458 dynParticle.SetKineticEnergy(kinEnergy);
459 currentModel->GetChargeSquareRatio(p, mat, kinEnergy);
460 currentModel->CorrectionsAlongStep(couple,&dynParticle,eloss,niel,length);
461 res = eloss/length;
462
463 if(verbose > 1) {
464 G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
465 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
466 << G4endl;
467 }
468 }
469 }
470 }
471 if(verbose > 0) {
472 G4cout << "Sum: E(MeV)= " << kinEnergy/MeV
473 << " DEDX(MeV/mm)= " << res*mm/MeV
474 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
475 << " cut(MeV)= " << cut/MeV
476 << " " << p->GetParticleName()
477 << " in " << currentMaterialName
478 << " Zi^2= " << chargeSquare
479 << " isIon=" << isIon
480 << G4endl;
481 }
482 }
483 return res;
484}
485
486//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487
488G4double G4EmCalculator::ComputeElectronicDEDX(G4double kinEnergy,
489 const G4ParticleDefinition* part,
490 const G4Material* mat,
491 G4double cut)
492{
493 SetupMaterial(mat);
494 G4double dedx = 0.0;
495 if(UpdateParticle(part, kinEnergy)) {
496
497 G4LossTableManager* lManager = G4LossTableManager::Instance();
498 const std::vector<G4VEnergyLossProcess*> vel =
499 lManager->GetEnergyLossProcessVector();
500 G4int n = vel.size();
501
502 //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
503 // << " n= " << n << G4endl;
504
505 for(G4int i=0; i<n; ++i) {
506 if(vel[i]) {
507 G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
508 if(ActiveForParticle(part, p)) {
509 //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
510 // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
511 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),mat,cut);
512 }
513 }
514 }
515 }
516 return dedx;
517}
518
519//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
520
521G4double G4EmCalculator::ComputeDEDXForCutInRange(G4double kinEnergy,
522 const G4ParticleDefinition* part,
523 const G4Material* mat,
524 G4double rangecut)
525{
526 SetupMaterial(mat);
527 G4double dedx = 0.0;
528 if(UpdateParticle(part, kinEnergy)) {
529
530 G4LossTableManager* lManager = G4LossTableManager::Instance();
531 const std::vector<G4VEnergyLossProcess*> vel =
532 lManager->GetEnergyLossProcessVector();
533 G4int n = vel.size();
534
535 if(mat != cutMaterial) {
536 cutMaterial = mat;
537 cutenergy[0] = ComputeEnergyCutFromRangeCut(rangecut, G4Gamma::Gamma(), mat);
538 cutenergy[1] = ComputeEnergyCutFromRangeCut(rangecut, G4Electron::Electron(), mat);
539 cutenergy[2] = ComputeEnergyCutFromRangeCut(rangecut, G4Positron::Positron(), mat);
540 }
541
542 //G4cout << "ComputeElectronicDEDX for " << part->GetParticleName()
543 // << " n= " << n << G4endl;
544
545 for(G4int i=0; i<n; ++i) {
546 if(vel[i]) {
547 G4VProcess* p = reinterpret_cast<G4VProcess*>(vel[i]);
548 if(ActiveForParticle(part, p)) {
549 //G4cout << "idx= " << i << " " << (vel[i])->GetProcessName()
550 // << " " << (vel[i])->Particle()->GetParticleName() << G4endl;
551 const G4ParticleDefinition* sec = (vel[i])->SecondaryParticle();
552 G4int idx = 0;
553 if(sec == G4Electron::Electron()) { idx = 1; }
554 else if(sec == G4Positron::Positron()) { idx = 2; }
555
556 dedx += ComputeDEDX(kinEnergy,part,(vel[i])->GetProcessName(),
557 mat,cutenergy[idx]);
558 }
559 }
560 }
561 }
562 return dedx;
563}
564
565//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
566
567G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy,
568 const G4ParticleDefinition* part,
569 const G4Material* mat,
570 G4double cut)
571{
572 G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
573 if(mass > 700.*MeV) { dedx += ComputeNuclearDEDX(kinEnergy,part,mat); }
574 return dedx;
575}
576
577//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
578
579G4double G4EmCalculator::ComputeNuclearDEDX(G4double kinEnergy,
580 const G4ParticleDefinition* p,
581 const G4Material* mat)
582{
583 G4double res = 0.0;
584 G4VEmProcess* nucst = FindDiscreteProcess(p, "nuclearStopping");
585 if(nucst) {
586 G4VEmModel* mod = nucst->GetModelByIndex();
587 if(mod) {
588 mod->SetFluctuationFlag(false);
589 res = mod->ComputeDEDXPerVolume(mat, p, kinEnergy);
590 }
591 }
592
593 if(verbose > 1) {
594 G4cout << p->GetParticleName() << " E(MeV)= " << kinEnergy/MeV
595 << " NuclearDEDX(MeV/mm)= " << res*mm/MeV
596 << " NuclearDEDX(MeV*cm^2/g)= "
597 << res*gram/(MeV*cm2*mat->GetDensity())
598 << G4endl;
599 }
600 return res;
601}
602
603//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
604
605G4double G4EmCalculator::ComputeCrossSectionPerVolume(
606 G4double kinEnergy,
607 const G4ParticleDefinition* p,
608 const G4String& processName,
609 const G4Material* mat,
610 G4double cut)
611{
612 SetupMaterial(mat);
613 G4double res = 0.0;
614 if(UpdateParticle(p, kinEnergy)) {
615 if(FindEmModel(p, processName, kinEnergy)) {
616 G4double e = kinEnergy;
617 G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
618 if(baseParticle) {
619 e *= kinEnergy*massRatio;
620 res = currentModel->CrossSectionPerVolume(
621 mat, baseParticle, e, aCut, e) * chargeSquare;
622 } else {
623 res = currentModel->CrossSectionPerVolume(mat, p, e, aCut, e);
624 }
625 if(verbose>0) {
626 G4cout << "G4EmCalculator::ComputeXSPerVolume: E(MeV)= " << kinEnergy/MeV
627 << " cross(cm-1)= " << res*cm
628 << " cut(keV)= " << aCut/keV
629 << " " << p->GetParticleName()
630 << " in " << mat->GetName()
631 << G4endl;
632 }
633 }
634 }
635 return res;
636}
637
638//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
639
640G4double G4EmCalculator::ComputeCrossSectionPerAtom(
641 G4double kinEnergy,
642 const G4ParticleDefinition* p,
643 const G4String& processName,
644 G4double Z, G4double A,
645 G4double cut)
646{
647 G4double res = 0.0;
648 if(UpdateParticle(p, kinEnergy)) {
649 G4int iz = G4lrint(Z);
650 CheckMaterial(iz);
651 if(FindEmModel(p, processName, kinEnergy)) {
652 G4double e = kinEnergy;
653 G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
654 if(baseParticle) {
655 e *= kinEnergy*massRatio;
656 currentModel->InitialiseForElement(baseParticle, iz);
657 res = currentModel->ComputeCrossSectionPerAtom(
658 baseParticle, e, Z, A, aCut) * chargeSquare;
659 } else {
660 currentModel->InitialiseForElement(p, iz);
661 res = currentModel->ComputeCrossSectionPerAtom(p, e, Z, A, aCut);
662 }
663 if(verbose>0) {
664 G4cout << "E(MeV)= " << kinEnergy/MeV
665 << " cross(barn)= " << res/barn
666 << " " << p->GetParticleName()
667 << " Z= " << Z << " A= " << A/(g/mole) << " g/mole"
668 << " cut(keV)= " << aCut/keV
669 << G4endl;
670 }
671 }
672 }
673 return res;
674}
675
676//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
677
678G4double G4EmCalculator::ComputeCrossSectionPerShell(G4double kinEnergy,
679 const G4ParticleDefinition* p,
680 const G4String& processName,
681 G4int Z, G4int shellIdx,
682 G4double cut)
683{
684 G4double res = 0.0;
685 if(UpdateParticle(p, kinEnergy)) {
686 CheckMaterial(Z);
687 if(FindEmModel(p, processName, kinEnergy)) {
688 G4double e = kinEnergy;
689 G4double aCut = std::max(cut, theParameters->LowestElectronEnergy());
690 if(baseParticle) {
691 e *= kinEnergy*massRatio;
692 currentModel->InitialiseForElement(baseParticle, Z);
693 res = currentModel->ComputeCrossSectionPerShell(baseParticle, Z, shellIdx,
694 e, aCut) * chargeSquare;
695 } else {
696 currentModel->InitialiseForElement(p, Z);
697 res = currentModel->ComputeCrossSectionPerAtom(p, Z, shellIdx, e, aCut);
698 }
699 if(verbose>0) {
700 G4cout << "E(MeV)= " << kinEnergy/MeV
701 << " cross(barn)= " << res/barn
702 << " " << p->GetParticleName()
703 << " Z= " << Z << " shellIdx= " << shellIdx
704 << " cut(keV)= " << aCut/keV
705 << G4endl;
706 }
707 }
708 }
709 return res;
710}
711
712//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
713
714G4double
715G4EmCalculator::ComputeGammaAttenuationLength(G4double kinEnergy,
716 const G4Material* mat)
717{
718 G4double res = 0.0;
719 const G4ParticleDefinition* gamma = G4Gamma::Gamma();
720 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "conv", mat, 0.0);
721 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "compt", mat, 0.0);
722 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "phot", mat, 0.0);
723 res += ComputeCrossSectionPerVolume(kinEnergy, gamma, "Rayl", mat, 0.0);
724 if(res > 0.0) { res = 1.0/res; }
725 return res;
726}
727
728//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
729
730G4double G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom(
731 const G4String& particle,
732 G4int Z,
733 G4AtomicShellEnumerator shell,
734 G4double kinEnergy,
735 const G4Material* mat)
736{
737 G4double res = 0.0;
738 const G4ParticleDefinition* p = FindParticle(particle);
739 G4VAtomDeexcitation* ad = manager->AtomDeexcitation();
740 if(p && ad) {
741 res = ad->ComputeShellIonisationCrossSectionPerAtom(p, Z, shell,
742 kinEnergy, mat);
743 }
744 return res;
745}
746
747//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
748
749G4double G4EmCalculator::ComputeMeanFreePath(G4double kinEnergy,
750 const G4ParticleDefinition* p,
751 const G4String& processName,
752 const G4Material* mat,
753 G4double cut)
754{
755 G4double mfp = DBL_MAX;
756 G4double x = ComputeCrossSectionPerVolume(kinEnergy, p, processName, mat, cut);
757 if(x > 0.0) { mfp = 1.0/x; }
758 if(verbose>1) {
759 G4cout << "E(MeV)= " << kinEnergy/MeV
760 << " MFP(mm)= " << mfp/mm
761 << " " << p->GetParticleName()
762 << " in " << mat->GetName()
763 << G4endl;
764 }
765 return mfp;
766}
767
768//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
769
770G4double G4EmCalculator::ComputeEnergyCutFromRangeCut(
771 G4double range,
772 const G4ParticleDefinition* part,
773 const G4Material* mat)
774{
775 return G4ProductionCutsTable::GetProductionCutsTable()->
776 ConvertRangeToEnergy(part, mat, range);
777}
778
779//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
780
781G4bool G4EmCalculator::UpdateParticle(const G4ParticleDefinition* p,
782 G4double kinEnergy)
783{
784 if(p != currentParticle) {
785
786 // new particle
787 currentParticle = p;
788 dynParticle.SetDefinition(const_cast<G4ParticleDefinition*>(p));
789 dynParticle.SetKineticEnergy(kinEnergy);
790 baseParticle = 0;
791 currentParticleName = p->GetParticleName();
792 massRatio = 1.0;
793 mass = p->GetPDGMass();
794 chargeSquare = 1.0;
795 currentProcess = FindEnergyLossProcess(p);
796 currentProcessName = "";
797 isIon = false;
798
799 // ionisation process exist
800 if(currentProcess) {
801 currentProcessName = currentProcess->GetProcessName();
802 baseParticle = currentProcess->BaseParticle();
803
804 // base particle is used
805 if(baseParticle) {
806 massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
807 G4double q = p->GetPDGCharge()/baseParticle->GetPDGCharge();
808 chargeSquare = q*q;
809 }
810
811 if(p->GetParticleType() == "nucleus"
812 && currentParticleName != "deuteron"
813 && currentParticleName != "triton"
814 && currentParticleName != "alpha+"
815 && currentParticleName != "helium"
816 && currentParticleName != "hydrogen"
817 ) {
818 isIon = true;
819 massRatio = theGenericIon->GetPDGMass()/p->GetPDGMass();
820 baseParticle = theGenericIon;
821 if(verbose>1) {
822 G4cout << "\n G4EmCalculator::UpdateParticle: isIon 1 "
823 << p->GetParticleName()
824 << " in " << currentMaterial->GetName()
825 << " e= " << kinEnergy << G4endl;
826 }
827 }
828 }
829 }
830
831 // Effective charge for ions
832 if(isIon) {
833 chargeSquare =
834 corr->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
835 * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
836 if(currentProcess) {
837 currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
838 if(verbose>1) {
839 G4cout <<"\n NewIon: massR= "<< massRatio << " q2= "
840 << chargeSquare << " " << currentProcess << G4endl;
841 }
842 }
843 }
844 return true;
845}
846
847//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
848
849const G4ParticleDefinition* G4EmCalculator::FindParticle(const G4String& name)
850{
851 const G4ParticleDefinition* p = 0;
852 if(name != currentParticleName) {
853 p = G4ParticleTable::GetParticleTable()->FindParticle(name);
854 if(!p) {
855 G4cout << "### WARNING: G4EmCalculator::FindParticle fails to find "
856 << name << G4endl;
857 }
858 } else {
859 p = currentParticle;
860 }
861 return p;
862}
863
864//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
865
866const G4ParticleDefinition* G4EmCalculator::FindIon(G4int Z, G4int A)
867{
868 const G4ParticleDefinition* p = ionTable->GetIon(Z,A,0);
869 return p;
870}
871
872//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
873
874const G4Material* G4EmCalculator::FindMaterial(const G4String& name)
875{
876 if(name != currentMaterialName) {
877 SetupMaterial(G4Material::GetMaterial(name, false));
878 if(!currentMaterial) {
879 G4cout << "### WARNING: G4EmCalculator::FindMaterial fails to find "
880 << name << G4endl;
881 }
882 }
883 return currentMaterial;
884}
885
886//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
887
888const G4Region* G4EmCalculator::FindRegion(const G4String& reg)
889{
890 const G4Region* r = 0;
891 if(reg != "" && reg != "world") {
892 r = G4RegionStore::GetInstance()->GetRegion(reg);
893 } else {
894 r = G4RegionStore::GetInstance()->GetRegion("DefaultRegionForTheWorld");
895 }
896 return r;
897}
898
899//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
900
901const G4MaterialCutsCouple* G4EmCalculator::FindCouple(
902 const G4Material* material,
903 const G4Region* region)
904{
905 const G4MaterialCutsCouple* couple = nullptr;
906 SetupMaterial(material);
907 if(currentMaterial) {
908 // Access to materials
909 const G4ProductionCutsTable* theCoupleTable=
910 G4ProductionCutsTable::GetProductionCutsTable();
911 const G4Region* r = region;
912 if(r) {
913 couple = theCoupleTable->GetMaterialCutsCouple(material,
914 r->GetProductionCuts());
915 } else {
916 G4RegionStore* store = G4RegionStore::GetInstance();
917 size_t nr = store->size();
918 if(0 < nr) {
919 for(size_t i=0; i<nr; ++i) {
920 couple = theCoupleTable->GetMaterialCutsCouple(
921 material, ((*store)[i])->GetProductionCuts());
922 if(couple) { break; }
923 }
924 }
925 }
926 }
927 if(!couple) {
928 G4ExceptionDescription ed;
929 ed << "G4EmCalculator::FindCouple: fail for material <"
930 << currentMaterialName << ">";
931 if(region) { ed << " and region " << region->GetName(); }
932 G4Exception("G4EmCalculator::FindCouple", "em0078",
933 FatalException, ed);
934 }
935 return couple;
936}
937
938//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
939
940G4bool G4EmCalculator::UpdateCouple(const G4Material* material, G4double cut)
941{
942 SetupMaterial(material);
943 if(!currentMaterial) { return false; }
944 for (G4int i=0; i<nLocalMaterials; ++i) {
945 if(material == localMaterials[i] && cut == localCuts[i]) {
946 currentCouple = localCouples[i];
947 currentCoupleIndex = currentCouple->GetIndex();
948 currentCut = cut;
949 return true;
950 }
951 }
952 const G4MaterialCutsCouple* cc = new G4MaterialCutsCouple(material);
953 localMaterials.push_back(material);
954 localCouples.push_back(cc);
955 localCuts.push_back(cut);
956 nLocalMaterials++;
957 currentCouple = cc;
958 currentCoupleIndex = currentCouple->GetIndex();
959 currentCut = cut;
960 return true;
961}
962
963//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
964
965void G4EmCalculator::FindLambdaTable(const G4ParticleDefinition* p,
966 const G4String& processName,
967 G4double kinEnergy, G4int& proctype)
968{
969 // Search for the process
970 if (!currentLambda || p != lambdaParticle || processName != lambdaName) {
971 lambdaName = processName;
972 currentLambda = nullptr;
973 lambdaParticle = p;
974
975 const G4ParticleDefinition* part = p;
976 if(isIon) { part = theGenericIon; }
977
978 // Search for energy loss process
979 currentName = processName;
980 currentModel = nullptr;
981 loweModel = nullptr;
982
983 G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
984 if(elproc) {
985 currentLambda = elproc->LambdaTable();
986 proctype = 0;
987 if(currentLambda) {
988 isApplicable = true;
989 if(verbose>1) {
990 G4cout << "G4VEnergyLossProcess is found out: " << currentName
991 << G4endl;
992 }
993 }
994 curProcess = elproc;
995 return;
996 }
997
998 // Search for discrete process
999 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1000 if(proc) {
1001 currentLambda = proc->LambdaTable();
1002 proctype = 1;
1003 if(currentLambda) {
1004 isApplicable = true;
1005 if(verbose>1) {
1006 G4cout << "G4VEmProcess is found out: " << currentName << G4endl;
1007 }
1008 }
1009 curProcess = proc;
1010 return;
1011 }
1012
1013 // Search for msc process
1014 G4VMultipleScattering* msc = FindMscProcess(part, processName);
1015 if(msc) {
1016 currentModel = msc->SelectModel(kinEnergy,0);
1017 proctype = 2;
1018 if(currentModel) {
1019 currentLambda = currentModel->GetCrossSectionTable();
1020 if(currentLambda) {
1021 isApplicable = true;
1022 if(verbose>1) {
1023 G4cout << "G4VMultipleScattering is found out: " << currentName
1024 << G4endl;
1025 }
1026 }
1027 }
1028 curProcess = msc;
1029 }
1030 }
1031}
1032
1033//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1034
1035G4bool G4EmCalculator::FindEmModel(const G4ParticleDefinition* p,
1036 const G4String& processName,
1037 G4double kinEnergy)
1038{
1039 isApplicable = false;
1040 if(!p || !currentMaterial) {
1041 G4cout << "G4EmCalculator::FindEmModel WARNING: no particle"
1042 << " or materail defined; particle: " << p << G4endl;
1043 return isApplicable;
1044 }
1045 G4String partname = p->GetParticleName();
1046 const G4ParticleDefinition* part = p;
1047 G4double scaledEnergy = kinEnergy*massRatio;
1048 if(isIon) { part = theGenericIon; }
1049
1050 if(verbose > 1) {
1051 G4cout << "## G4EmCalculator::FindEmModel for " << partname
1052 << " (type= " << p->GetParticleType()
1053 << ") and " << processName << " at E(MeV)= " << scaledEnergy
1054 << G4endl;
1055 if(p != part) { G4cout << " GenericIon is the base particle" << G4endl; }
1056 }
1057
1058 // Search for energy loss process
1059 currentName = processName;
1060 currentModel = nullptr;
1061 loweModel = nullptr;
1062 size_t idx = 0;
1063
1064 G4VEnergyLossProcess* elproc = FindEnLossProcess(part, processName);
1065 if(elproc) {
1066 currentModel = elproc->SelectModelForMaterial(scaledEnergy, idx);
1067 currentModel->InitialiseForMaterial(part, currentMaterial);
1068 currentModel->SetupForMaterial(part, currentMaterial, scaledEnergy);
1069 G4double eth = currentModel->LowEnergyLimit();
1070 if(eth > 0.0) {
1071 loweModel = elproc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1072 if(loweModel == currentModel) { loweModel = nullptr; }
1073 else {
1074 loweModel->InitialiseForMaterial(part, currentMaterial);
1075 loweModel->SetupForMaterial(part, currentMaterial, eth - CLHEP::eV);
1076 }
1077 }
1078 }
1079
1080 // Search for discrete process
1081 if(!currentModel) {
1082 G4VEmProcess* proc = FindDiscreteProcess(part, processName);
1083 if(proc) {
1084 currentModel = proc->SelectModelForMaterial(kinEnergy, idx);
1085 currentModel->InitialiseForMaterial(part, currentMaterial);
1086 currentModel->SetupForMaterial(part, currentMaterial, kinEnergy);
1087 G4double eth = currentModel->LowEnergyLimit();
1088 if(eth > 0.0) {
1089 loweModel = proc->SelectModelForMaterial(eth - CLHEP::eV, idx);
1090 if(loweModel == currentModel) { loweModel = nullptr; }
1091 else {
1092 loweModel->InitialiseForMaterial(part, currentMaterial);
1093 loweModel->SetupForMaterial(part, currentMaterial, eth - CLHEP::eV);
1094 }
1095 }
1096 }
1097 }
1098
1099 // Search for msc process
1100 if(!currentModel) {
1101 G4VMultipleScattering* proc = FindMscProcess(part, processName);
1102 if(proc) {
1103 currentModel = proc->SelectModel(kinEnergy, idx);
1104 loweModel = nullptr;
1105 }
1106 }
1107 if(currentModel) {
1108 if(loweModel == currentModel) { loweModel = nullptr; }
1109 isApplicable = true;
1110 currentModel->InitialiseForMaterial(part, currentMaterial);
1111 if(loweModel) {
1112 loweModel->InitialiseForMaterial(part, currentMaterial);
1113 }
1114 if(verbose > 1) {
1115 G4cout << " Model <" << currentModel->GetName()
1116 << "> Emin(MeV)= " << currentModel->LowEnergyLimit()/MeV
1117 << " for " << part->GetParticleName();
1118 if(elproc) {
1119 G4cout << " and " << elproc->GetProcessName() << " " << elproc
1120 << G4endl;
1121 }
1122 if(loweModel) {
1123 G4cout << " LowEnergy model <" << loweModel->GetName() << ">";
1124 }
1125 G4cout << G4endl;
1126 }
1127 }
1128 return isApplicable;
1129}
1130
1131//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1132
1133G4VEnergyLossProcess* G4EmCalculator::FindEnergyLossProcess(
1134 const G4ParticleDefinition* p)
1135{
1136 G4VEnergyLossProcess* elp = nullptr;
1137 G4String partname = p->GetParticleName();
1138 const G4ParticleDefinition* part = p;
1139
1140 if(p->GetParticleType() == "nucleus"
1141 && currentParticleName != "deuteron"
1142 && currentParticleName != "triton"
1143 && currentParticleName != "He3"
1144 && currentParticleName != "alpha"
1145 && currentParticleName != "alpha+"
1146 && currentParticleName != "helium"
1147 && currentParticleName != "hydrogen"
1148 ) { part = theGenericIon; }
1149
1150 elp = manager->GetEnergyLossProcess(part);
1151 /*
1152 G4cout << "\n G4EmCalculator::FindEnergyLossProcess: for " << p->GetParticleName()
1153 << " found " << elp->GetProcessName() << " of "
1154 << elp->Particle()->GetParticleName() << " " << elp << G4endl;
1155 */
1156 return elp;
1157}
1158
1159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1160
1161G4VEnergyLossProcess*
1162G4EmCalculator::FindEnLossProcess(const G4ParticleDefinition* part,
1163 const G4String& processName)
1164{
1165 G4VEnergyLossProcess* proc = 0;
1166 const std::vector<G4VEnergyLossProcess*> v =
1167 manager->GetEnergyLossProcessVector();
1168 G4int n = v.size();
1169 for(G4int i=0; i<n; ++i) {
1170 if((v[i])->GetProcessName() == processName) {
1171 G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1172 if(ActiveForParticle(part, p)) {
1173 proc = v[i];
1174 break;
1175 }
1176 }
1177 }
1178 return proc;
1179}
1180
1181//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1182
1183G4VEmProcess*
1184G4EmCalculator::FindDiscreteProcess(const G4ParticleDefinition* part,
1185 const G4String& processName)
1186{
1187 G4VEmProcess* proc = nullptr;
1188 auto v = manager->GetEmProcessVector();
1189 G4int n = v.size();
1190 for(G4int i=0; i<n; ++i) {
1191 auto pName = v[i]->GetProcessName();
1192 if(pName == "GammaGeneralProc") {
1193 proc = v[i]->GetEmProcess(processName);
1194 break;
1195 } else if(pName == processName) {
1196 G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1197 if(ActiveForParticle(part, p)) {
1198 proc = v[i];
1199 break;
1200 }
1201 }
1202 }
1203 return proc;
1204}
1205
1206//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1207
1208G4VMultipleScattering*
1209G4EmCalculator::FindMscProcess(const G4ParticleDefinition* part,
1210 const G4String& processName)
1211{
1212 G4VMultipleScattering* proc = 0;
1213 const std::vector<G4VMultipleScattering*> v =
1214 manager->GetMultipleScatteringVector();
1215 G4int n = v.size();
1216 for(G4int i=0; i<n; ++i) {
1217 if((v[i])->GetProcessName() == processName) {
1218 G4VProcess* p = reinterpret_cast<G4VProcess*>(v[i]);
1219 if(ActiveForParticle(part, p)) {
1220 proc = v[i];
1221 break;
1222 }
1223 }
1224 }
1225 return proc;
1226}
1227
1228//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1229
1230G4VProcess* G4EmCalculator::FindProcess(const G4ParticleDefinition* part,
1231 const G4String& processName)
1232{
1233 G4VProcess* proc = 0;
1234 const G4ProcessManager* procman = part->GetProcessManager();
1235 G4ProcessVector* pv = procman->GetProcessList();
1236 G4int nproc = pv->size();
1237 for(G4int i=0; i<nproc; ++i) {
1238 if(processName == (*pv)[i]->GetProcessName()) {
1239 proc = (*pv)[i];
1240 break;
1241 }
1242 }
1243 return proc;
1244}
1245
1246//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1247
1248G4bool G4EmCalculator::ActiveForParticle(const G4ParticleDefinition* part,
1249 G4VProcess* proc)
1250{
1251 G4ProcessManager* pm = part->GetProcessManager();
1252 G4ProcessVector* pv = pm->GetProcessList();
1253 G4int n = pv->size();
1254 G4bool res = false;
1255 for(G4int i=0; i<n; ++i) {
1256 if((*pv)[i] == proc) {
1257 if(pm->GetProcessActivation(i)) { res = true; }
1258 break;
1259 }
1260 }
1261 return res;
1262}
1263
1264//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1265
1266void G4EmCalculator::SetupMaterial(const G4Material* mat)
1267{
1268 if(mat) {
1269 currentMaterial = mat;
1270 currentMaterialName = mat->GetName();
1271 } else {
1272 currentMaterial = 0;
1273 currentMaterialName = "";
1274 }
1275}
1276
1277//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1278
1279void G4EmCalculator::SetupMaterial(const G4String& mname)
1280{
1281 SetupMaterial(nist->FindOrBuildMaterial(mname));
1282}
1283
1284//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1285
1286void G4EmCalculator::CheckMaterial(G4int Z)
1287{
1288 G4bool isFound = false;
1289 if(currentMaterial) {
1290 size_t nn = currentMaterial->GetNumberOfElements();
1291 for(size_t i=0; i<nn; ++i) {
1292 if(Z == currentMaterial->GetElement(i)->GetZasInt()) {
1293 isFound = true;
1294 break;
1295 }
1296 }
1297 }
1298 if(!isFound) {
1299 SetupMaterial(nist->FindOrBuildSimpleMaterial(Z));
1300 }
1301}
1302
1303//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1304
1305void G4EmCalculator::SetVerbose(G4int verb)
1306{
1307 verbose = verb;
1308}
1309
1310//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1311
G4AtomicShellEnumerator
G4AtomicShellEnumerator
Definition: G4AtomicShellEnumerator.hh:50
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:60
JustWarning
@ JustWarning
Definition: G4ExceptionSeverity.hh:65
FatalException
@ FatalException
Definition: G4ExceptionSeverity.hh:61
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4double
double G4double
Definition: G4Types.hh:83
G4bool
bool G4bool
Definition: G4Types.hh:86
G4int
int G4int
Definition: G4Types.hh:85
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout
G4DynamicParticle::SetDefinition
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
Definition: G4DynamicParticle.cc:325
G4DynamicParticle::SetKineticEnergy
void SetKineticEnergy(G4double aEnergy)
G4Electron::Electron
static G4Electron * Electron()
Definition: G4Electron.cc:93
G4Element::GetZasInt
G4int GetZasInt() const
Definition: G4Element.hh:131
G4EmCalculator::FindParticle
const G4ParticleDefinition * FindParticle(const G4String &)
Definition: G4EmCalculator.cc:849
G4EmCalculator::ComputeMeanFreePath
G4double ComputeMeanFreePath(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=0.0)
Definition: G4EmCalculator.cc:749
G4EmCalculator::GetShellIonisationCrossSectionPerAtom
G4double GetShellIonisationCrossSectionPerAtom(const G4String &part, G4int Z, G4AtomicShellEnumerator shell, G4double kinEnergy)
Definition: G4EmCalculator.cc:306
G4EmCalculator::GetDEDX
G4double GetDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:121
G4EmCalculator::ComputeTotalDEDX
G4double ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, G4double cut=DBL_MAX)
Definition: G4EmCalculator.cc:567
G4EmCalculator::GetCSDARange
G4double GetCSDARange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:185
G4EmCalculator::GetRange
G4double GetRange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:216
G4EmCalculator::ComputeNuclearDEDX
G4double ComputeNuclearDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *)
Definition: G4EmCalculator.cc:579
G4EmCalculator::FindIon
const G4ParticleDefinition * FindIon(G4int Z, G4int A)
Definition: G4EmCalculator.cc:866
G4EmCalculator::ComputeGammaAttenuationLength
G4double ComputeGammaAttenuationLength(G4double kinEnergy, const G4Material *)
Definition: G4EmCalculator.cc:715
G4EmCalculator::ComputeDEDX
G4double ComputeDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=DBL_MAX)
Definition: G4EmCalculator.cc:372
G4EmCalculator::SetVerbose
void SetVerbose(G4int val)
Definition: G4EmCalculator.cc:1305
G4EmCalculator::ComputeEnergyCutFromRangeCut
G4double ComputeEnergyCutFromRangeCut(G4double range, const G4ParticleDefinition *, const G4Material *)
Definition: G4EmCalculator.cc:770
G4EmCalculator::GetKinEnergy
G4double GetKinEnergy(G4double range, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:232
G4EmCalculator::GetRangeFromRestricteDEDX
G4double GetRangeFromRestricteDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:162
G4EmCalculator::ComputeCrossSectionPerVolume
G4double ComputeCrossSectionPerVolume(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, G4double cut=0.0)
Definition: G4EmCalculator.cc:605
G4EmCalculator::ComputeCrossSectionPerAtom
G4double ComputeCrossSectionPerAtom(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, G4double Z, G4double A, G4double cut=0.0)
Definition: G4EmCalculator.cc:640
G4EmCalculator::PrintDEDXTable
void PrintDEDXTable(const G4ParticleDefinition *)
Definition: G4EmCalculator.cc:344
G4EmCalculator::FindRegion
const G4Region * FindRegion(const G4String &)
Definition: G4EmCalculator.cc:888
G4EmCalculator::FindProcess
G4VProcess * FindProcess(const G4ParticleDefinition *part, const G4String &processName)
Definition: G4EmCalculator.cc:1230
G4EmCalculator::GetCrossSectionPerVolume
G4double GetCrossSectionPerVolume(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:254
G4EmCalculator::ComputeCrossSectionPerShell
G4double ComputeCrossSectionPerShell(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, G4int Z, G4int shellIdx, G4double cut=0.0)
Definition: G4EmCalculator.cc:678
G4EmCalculator::ComputeShellIonisationCrossSectionPerAtom
G4double ComputeShellIonisationCrossSectionPerAtom(const G4String &part, G4int Z, G4AtomicShellEnumerator shell, G4double kinEnergy, const G4Material *mat=nullptr)
Definition: G4EmCalculator.cc:730
G4EmCalculator::FindCouple
const G4MaterialCutsCouple * FindCouple(const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:901
G4EmCalculator::SetupMaterial
void SetupMaterial(const G4Material *)
Definition: G4EmCalculator.cc:1266
G4EmCalculator::PrintRangeTable
void PrintRangeTable(const G4ParticleDefinition *)
Definition: G4EmCalculator.cc:353
G4EmCalculator::PrintInverseRangeTable
void PrintInverseRangeTable(const G4ParticleDefinition *)
Definition: G4EmCalculator.cc:362
G4EmCalculator::ComputeDEDXForCutInRange
G4double ComputeDEDXForCutInRange(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *mat, G4double rangecut=DBL_MAX)
Definition: G4EmCalculator.cc:521
G4EmCalculator::GetMeanFreePath
G4double GetMeanFreePath(G4double kinEnergy, const G4ParticleDefinition *, const G4String &processName, const G4Material *, const G4Region *r=nullptr)
Definition: G4EmCalculator.cc:323
G4EmCalculator::ComputeElectronicDEDX
G4double ComputeElectronicDEDX(G4double kinEnergy, const G4ParticleDefinition *, const G4Material *mat, G4double cut=DBL_MAX)
Definition: G4EmCalculator.cc:488
G4EmCalculator::FindMaterial
const G4Material * FindMaterial(const G4String &)
Definition: G4EmCalculator.cc:874
G4EmCorrections::EffectiveChargeSquareRatio
G4double EffectiveChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.hh:317
G4EmCorrections::EffectiveChargeCorrection
G4double EffectiveChargeCorrection(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4EmCorrections.cc:766
G4EmParameters::Instance
static G4EmParameters * Instance()
Definition: G4EmParameters.cc:67
G4EmParameters::BuildCSDARange
G4bool BuildCSDARange() const
Definition: G4EmParameters.cc:193
G4EmParameters::LowestElectronEnergy
G4double LowestElectronEnergy() const
Definition: G4EmParameters.cc:600
G4Gamma::Gamma
static G4Gamma * Gamma()
Definition: G4Gamma.cc:85
G4GenericIon::GenericIon
static G4GenericIon * GenericIon()
Definition: G4GenericIon.cc:92
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:87
G4LossTableManager::GetEmProcessVector
const std::vector< G4VEmProcess * > & GetEmProcessVector()
Definition: G4LossTableManager.cc:994
G4LossTableManager::GetCSDARange
G4double GetCSDARange(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:333
G4LossTableManager::GetEnergyLossProcess
G4VEnergyLossProcess * GetEnergyLossProcess(const G4ParticleDefinition *)
Definition: G4LossTableManager.cc:416
G4LossTableManager::GetEnergy
G4double GetEnergy(const G4ParticleDefinition *aParticle, G4double range, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:367
G4LossTableManager::GetMultipleScatteringVector
const std::vector< G4VMultipleScattering * > & GetMultipleScatteringVector()
Definition: G4LossTableManager.cc:1002
G4LossTableManager::GetEnergyLossProcessVector
const std::vector< G4VEnergyLossProcess * > & GetEnergyLossProcessVector()
Definition: G4LossTableManager.cc:987
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.hh:404
G4LossTableManager::GetRangeFromRestricteDEDX
G4double GetRangeFromRestricteDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:344
G4LossTableManager::EmCorrections
G4EmCorrections * EmCorrections()
Definition: G4LossTableManager.hh:397
G4LossTableManager::GetDEDX
G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double kineticEnergy, const G4MaterialCutsCouple *couple)
Definition: G4LossTableManager.hh:311
G4Material::GetDensity
G4double GetDensity() const
Definition: G4Material.hh:178
G4Material::GetElement
const G4Element * GetElement(G4int iel) const
Definition: G4Material.hh:200
G4Material::GetNumberOfElements
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
G4Material::GetName
const G4String & GetName() const
Definition: G4Material.hh:175
G4Material::GetMaterial
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:651
G4NistManager::FindOrBuildSimpleMaterial
G4Material * FindOrBuildSimpleMaterial(G4int Z, G4bool warning=false)
Definition: G4NistManager.hh:509
G4NistManager::FindOrBuildMaterial
G4Material * FindOrBuildMaterial(const G4String &name, G4bool isotopes=true, G4bool warning=false)
Definition: G4NistManager.hh:491
G4NistManager::Instance
static G4NistManager * Instance()
Definition: G4NistManager.cc:70
G4ParticleDefinition::GetProcessManager
G4ProcessManager * GetProcessManager() const
Definition: G4ParticleDefinition.cc:237
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:129
G4ParticleDefinition::GetPDGCharge
G4double GetPDGCharge() const
Definition: G4ParticleDefinition.hh:112
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:108
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:695
G4ParticleTable::FindParticle
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
Definition: G4ParticleTable.cc:586
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:87
G4Positron::Positron
static G4Positron * Positron()
Definition: G4Positron.cc:93
G4ProcessManager::GetProcessActivation
G4bool GetProcessActivation(G4VProcess *aProcess) const
Definition: G4ProcessManager.cc:1210
G4ProcessManager::GetProcessList
G4ProcessVector * GetProcessList() const
G4ProcessVector::size
std::size_t size() const
G4ProductionCutsTable::GetMaterialCutsCouple
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
Definition: G4ProductionCutsTable.hh:246
G4ProductionCutsTable::GetProductionCutsTable
static G4ProductionCutsTable * GetProductionCutsTable()
Definition: G4ProductionCutsTable.cc:58
G4RegionStore::GetInstance
static G4RegionStore * GetInstance()
Definition: G4RegionStore.cc:156
G4RegionStore::GetRegion
G4Region * GetRegion(const G4String &name, G4bool verbose=true) const
Definition: G4RegionStore.cc:211
G4Region::GetProductionCuts
G4ProductionCuts * GetProductionCuts() const
G4Region::GetName
const G4String & GetName() const
G4VAtomDeexcitation::ComputeShellIonisationCrossSectionPerAtom
virtual G4double ComputeShellIonisationCrossSectionPerAtom(const G4ParticleDefinition *, G4int Z, G4AtomicShellEnumerator shell, G4double kinE, const G4Material *mat=nullptr)=0
G4VAtomDeexcitation::GetShellIonisationCrossSectionPerAtom
virtual G4double GetShellIonisationCrossSectionPerAtom(const G4ParticleDefinition *, G4int Z, G4AtomicShellEnumerator shell, G4double kinE, const G4Material *mat=nullptr)=0
G4VEmModel::ComputeCrossSectionPerAtom
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:359
G4VEmModel::CorrectionsAlongStep
virtual void CorrectionsAlongStep(const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
Definition: G4VEmModel.cc:408
G4VEmModel::InitialiseForMaterial
virtual void InitialiseForMaterial(const G4ParticleDefinition *, const G4Material *)
Definition: G4VEmModel.cc:226
G4VEmModel::LowEnergyLimit
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:652
G4VEmModel::CrossSectionPerVolume
virtual G4double CrossSectionPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:254
G4VEmModel::SetCurrentCouple
void SetCurrentCouple(const G4MaterialCutsCouple *)
Definition: G4VEmModel.hh:465
G4VEmModel::SetFluctuationFlag
void SetFluctuationFlag(G4bool val)
Definition: G4VEmModel.hh:722
G4VEmModel::InitialiseForElement
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
Definition: G4VEmModel.cc:240
G4VEmModel::SetupForMaterial
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:449
G4VEmModel::ComputeCrossSectionPerShell
virtual G4double ComputeCrossSectionPerShell(const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.cc:369
G4VEmModel::GetChargeSquareRatio
virtual G4double GetChargeSquareRatio(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:391
G4VEmModel::GetName
const G4String & GetName() const
Definition: G4VEmModel.hh:827
G4VEmModel::ComputeDEDXPerVolume
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
Definition: G4VEmModel.cc:245
G4VEmModel::GetCrossSectionTable
G4PhysicsTable * GetCrossSectionTable()
Definition: G4VEmModel.hh:860
G4VEmProcess::GetModelByIndex
G4VEmModel * GetModelByIndex(G4int idx=0, G4bool ver=false) const
Definition: G4VEmProcess.cc:247
G4VEmProcess::CrossSectionPerVolume
G4double CrossSectionPerVolume(G4double kineticEnergy, const G4MaterialCutsCouple *couple, G4double logKinEnergy=DBL_MAX)
Definition: G4VEmProcess.cc:1014
G4VEmProcess::LambdaTable
G4PhysicsTable * LambdaTable() const
Definition: G4VEmProcess.hh:605
G4VEmProcess::SelectModelForMaterial
G4VEmModel * SelectModelForMaterial(G4double kinEnergy, size_t idxRegion) const
Definition: G4VEmProcess.hh:502
G4VEmProcess::GetEmProcess
virtual G4VEmProcess * GetEmProcess(const G4String &name)
Definition: G4VEmProcess.cc:1251
G4VEnergyLossProcess::BaseParticle
const G4ParticleDefinition * BaseParticle() const
Definition: G4VEnergyLossProcess.hh:980
G4VEnergyLossProcess::RangeTableForLoss
G4PhysicsTable * RangeTableForLoss() const
Definition: G4VEnergyLossProcess.hh:1109
G4VEnergyLossProcess::InverseRangeTable
G4PhysicsTable * InverseRangeTable() const
Definition: G4VEnergyLossProcess.hh:1116
G4VEnergyLossProcess::SelectModelForMaterial
G4VEmModel * SelectModelForMaterial(G4double kinEnergy, size_t &idx) const
Definition: G4VEnergyLossProcess.hh:610
G4VEnergyLossProcess::SetDynamicMassCharge
void SetDynamicMassCharge(G4double massratio, G4double charge2ratio)
Definition: G4VEnergyLossProcess.hh:635
G4VEnergyLossProcess::LambdaTable
G4PhysicsTable * LambdaTable() const
Definition: G4VEnergyLossProcess.hh:1123
G4VEnergyLossProcess::DEDXTable
G4PhysicsTable * DEDXTable() const
Definition: G4VEnergyLossProcess.hh:1060
G4VMscModel::GetTransportMeanFreePath
G4double GetTransportMeanFreePath(const G4ParticleDefinition *part, G4double kinEnergy)
Definition: G4VMscModel.hh:405
G4VMultipleScattering::SelectModel
G4VEmModel * SelectModel(G4double kinEnergy, size_t idx)
Definition: G4VMultipleScattering.hh:310
G4VProcess::GetProcessName
const G4String & GetProcessName() const
Definition: G4VProcess.hh:382
G4lrint
int G4lrint(double ad)
Definition: templates.hh:134
DBL_MAX
#define DBL_MAX
Definition: templates.hh:62