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