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