Geant4  10.03
G4PenelopeIonisationModel.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: G4PenelopeIonisationModel.cc 99415 2016-09-21 09:05:43Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 //
30 // History:
31 // --------
32 // 27 Jul 2010 L Pandola First complete implementation
33 // 18 Jan 2011 L.Pandola Stricter check on production of sub-treshold delta-rays.
34 // Should never happen now
35 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is active.
36 // Make sure that fluorescence/Auger is generated only if
37 // above threshold
38 // 25 May 2011 L Pandola Renamed (make v2008 as default Penelope)
39 // 26 Jan 2012 L Pandola Migration of AtomicDeexcitation to the new interface
40 // 09 Mar 2012 L Pandola Moved the management and calculation of
41 // cross sections to a separate class. Use a different method to
42 // get normalized shell cross sections
43 // 07 Oct 2013 L. Pandola Migration to MT
44 // 23 Jun 2015 L. Pandola Keep track of the PIXE flag, to avoid double-production of
45 // atomic de-excitation (bug #1761)
46 //
47 
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4ParticleDefinition.hh"
52 #include "G4MaterialCutsCouple.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4DynamicParticle.hh"
56 #include "G4AtomicShell.hh"
57 #include "G4Gamma.hh"
58 #include "G4Electron.hh"
59 #include "G4Positron.hh"
61 #include "G4PenelopeOscillator.hh"
63 #include "G4PhysicsFreeVector.hh"
64 #include "G4PhysicsLogVector.hh"
65 #include "G4LossTableManager.hh"
67 #include "G4EmParameters.hh"
68 #include "G4AutoLock.hh"
69 
70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71 namespace { G4Mutex PenelopeIonisationModelMutex = G4MUTEX_INITIALIZER; }
72 
74  const G4String& nam)
75  :G4VEmModel(nam),fParticleChange(0),fParticle(0),isInitialised(false),
76  fAtomDeexcitation(0),fPIXEflag(false),kineticEnergy1(0.*eV),
77  cosThetaPrimary(1.0),energySecondary(0.*eV),
78  cosThetaSecondary(0.0),targetOscillator(-1),
79  theCrossSectionHandler(0),fLocalTable(false)
80 {
83  // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
85  nBins = 200;
86 
87  if (part)
88  SetParticle(part);
89 
90  //
92  //
93  verboseLevel= 0;
94 
95  // Verbosity scale:
96  // 0 = nothing
97  // 1 = warning for energy non-conservation
98  // 2 = details of energy budget
99  // 3 = calculation of cross sections, file openings, sampling of atoms
100  // 4 = entering in methods
101 
102  // Atomic deexcitation model activated by default
103  SetDeexcitationFlag(true);
104 }
105 
106 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107 
109 {
110  if (IsMaster() || fLocalTable)
111  {
113  delete theCrossSectionHandler;
114  }
115 }
116 
117 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118 
120  const G4DataVector& theCuts)
121 {
122  if (verboseLevel > 3)
123  G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
124 
126  //Issue warning if the AtomicDeexcitation has not been declared
127  if (!fAtomDeexcitation)
128  {
129  G4cout << G4endl;
130  G4cout << "WARNING from G4PenelopeIonisationModel " << G4endl;
131  G4cout << "Atomic de-excitation module is not instantiated, so there will not be ";
132  G4cout << "any fluorescence/Auger emission." << G4endl;
133  G4cout << "Please make sure this is intended" << G4endl;
134  }
135 
136  if (fAtomDeexcitation)
138 
139  //If the PIXE flag is active, the PIXE interface will take care of the
140  //atomic de-excitation. The model does not need to do that.
141  //Issue warnings here
142  if (fPIXEflag && IsMaster() && particle==G4Electron::Electron())
143  {
145  G4cout << "======================================================================" << G4endl;
146  G4cout << "The G4PenelopeIonisationModel is being used with the PIXE flag ON." << G4endl;
147  G4cout << "Atomic de-excitation will be produced statistically by the PIXE " << G4endl;
148  G4cout << "interface by using the shell cross section --> " << theModel << G4endl;
149  G4cout << "The built-in model procedure for atomic de-excitation is disabled. " << G4endl;
150  G4cout << "*Please be sure this is intended*, or disable PIXE by" << G4endl;
151  G4cout << "/process/em/pixe false" << G4endl;
152  /*
153  if (theModel != "Penelope")
154  {
155  ed << "To use the PIXE interface with the Penelope-based shell cross sections" << G4endl;
156  ed << "/process/em/pixeElecXSmodel Penelope" << G4endl;
157 
158  }
159  */
160  G4cout << "======================================================================" << G4endl;
161 
162  }
163 
164 
165 
166  SetParticle(particle);
167 
168  //Only the master model creates/manages the tables. All workers get the
169  //pointer to the table, and use it as readonly
170  if (IsMaster() && particle == fParticle)
171  {
172 
173  //Set the number of bins for the tables. 20 points per decade
174  nBins = (size_t) (20*std::log10(HighEnergyLimit()/LowEnergyLimit()));
175  nBins = std::max(nBins,(size_t)100);
176 
177  //Clear and re-build the tables
179  {
180  delete theCrossSectionHandler;
182  }
185 
186  //Build tables for all materials
187  G4ProductionCutsTable* theCoupleTable =
189  for (size_t i=0;i<theCoupleTable->GetTableSize();i++)
190  {
191  const G4Material* theMat =
192  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
193  //Forces the building of the cross section tables
194  theCrossSectionHandler->BuildXSTable(theMat,theCuts.at(i),particle,
195  IsMaster());
196 
197  }
198 
199  if (verboseLevel > 2) {
200  G4cout << "Penelope Ionisation model v2008 is initialized " << G4endl
201  << "Energy range: "
202  << LowEnergyLimit() / keV << " keV - "
203  << HighEnergyLimit() / GeV << " GeV. Using "
204  << nBins << " bins."
205  << G4endl;
206  }
207  }
208 
209  if(isInitialised)
210  return;
212  isInitialised = true;
213 
214  return;
215 }
216 
217 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
218 
220  G4VEmModel *masterModel)
221 {
222  if (verboseLevel > 3)
223  G4cout << "Calling G4PenelopeIonisationModel::InitialiseLocal()" << G4endl;
224 
225  //
226  //Check that particle matches: one might have multiple master models (e.g.
227  //for e+ and e-).
228  //
229  if (part == fParticle)
230  {
231  //Get the const table pointers from the master to the workers
232  const G4PenelopeIonisationModel* theModel =
233  static_cast<G4PenelopeIonisationModel*> (masterModel);
234 
235  //Copy pointers to the data tables
237 
238  //copy data
239  nBins = theModel->nBins;
240 
241  //Same verbosity for all workers, as the master
242  verboseLevel = theModel->verboseLevel;
243  }
244 
245  return;
246 }
247 
248 
249 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250 
252  const G4ParticleDefinition*
253  theParticle,
255  G4double cutEnergy,
256  G4double)
257 {
258  // Penelope model v2008 to calculate the cross section for inelastic collisions above the
259  // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from
260  // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
261  //
262  // The total cross section is calculated analytically by taking
263  // into account the atomic oscillators coming into the play for a given threshold.
264  //
265  // For incident e- the maximum energy allowed for the delta-rays is energy/2.
266  // because particles are undistinghishable.
267  //
268  // The contribution is splitted in three parts: distant longitudinal collisions,
269  // distant transverse collisions and close collisions. Each term is described by
270  // its own analytical function.
271  // Fermi density correction is calculated analytically according to
272  // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
273  //
274  if (verboseLevel > 3)
275  G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl;
276 
277  SetupForMaterial(theParticle, material, energy);
278 
279  G4double totalCross = 0.0;
280  G4double crossPerMolecule = 0.;
281 
282  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
283  //not invoked
285  {
286  //create a **thread-local** version of the table. Used only for G4EmCalculator and
287  //Unit Tests
288  fLocalTable = true;
290  }
291 
292  const G4PenelopeCrossSection* theXS =
294  material,
295  cutEnergy);
296  if (!theXS)
297  {
298  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
299  //not filled up. This can happen in a UnitTest or via G4EmCalculator
300  if (verboseLevel > 0)
301  {
302  //Issue a G4Exception (warning) only in verbose mode
304  ed << "Unable to retrieve the cross section table for " << theParticle->GetParticleName() <<
305  " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl;
306  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
307  G4Exception("G4PenelopeIonisationModel::CrossSectionPerVolume()",
308  "em2038",JustWarning,ed);
309  }
310  //protect file reading via autolock
311  G4AutoLock lock(&PenelopeIonisationModelMutex);
312  theCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle);
313  lock.unlock();
314  //now it should be ok
315  theXS =
317  material,
318  cutEnergy);
319  }
320 
321 
322  if (theXS)
323  crossPerMolecule = theXS->GetHardCrossSection(energy);
324 
325  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
326  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
327 
328  if (verboseLevel > 3)
329  G4cout << "Material " << material->GetName() << " has " << atPerMol <<
330  "atoms per molecule" << G4endl;
331 
332  G4double moleculeDensity = 0.;
333  if (atPerMol)
334  moleculeDensity = atomDensity/atPerMol;
335 
336  G4double crossPerVolume = crossPerMolecule*moleculeDensity;
337 
338  if (verboseLevel > 2)
339  {
340  G4cout << "G4PenelopeIonisationModel " << G4endl;
341  G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " <<
342  energy/keV << " keV = " << (1./crossPerVolume)/mm << " mm" << G4endl;
343  if (theXS)
344  totalCross = (theXS->GetTotalCrossSection(energy))*moleculeDensity;
345  G4cout << "Total free path for ionisation (no threshold) at " <<
346  energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
347  }
348  return crossPerVolume;
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 
353 //This is a dummy method. Never inkoved by the tracking, it just issues
354 //a warning if one tries to get Cross Sections per Atom via the
355 //G4EmCalculator.
357  G4double,
358  G4double,
359  G4double,
360  G4double,
361  G4double)
362 {
363  G4cout << "*** G4PenelopeIonisationModel -- WARNING ***" << G4endl;
364  G4cout << "Penelope Ionisation model v2008 does not calculate cross section _per atom_ " << G4endl;
365  G4cout << "so the result is always zero. For physics values, please invoke " << G4endl;
366  G4cout << "GetCrossSectionPerVolume() or GetMeanFreePath() via the G4EmCalculator" << G4endl;
367  return 0;
368 }
369 
370 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
371 
373  const G4ParticleDefinition* theParticle,
374  G4double kineticEnergy,
375  G4double cutEnergy)
376 {
377  // Penelope model v2008 to calculate the stopping power for soft inelastic collisions
378  // below the threshold. It makes use of the Generalised Oscillator Strength (GOS)
379  // model from
380  // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
381  //
382  // The stopping power is calculated analytically using the dsigma/dW cross
383  // section from the GOS models, which includes separate contributions from
384  // distant longitudinal collisions, distant transverse collisions and
385  // close collisions. Only the atomic oscillators that come in the play
386  // (according to the threshold) are considered for the calculation.
387  // Differential cross sections have a different form for e+ and e-.
388  //
389  // Fermi density correction is calculated analytically according to
390  // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
391 
392  if (verboseLevel > 3)
393  G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl;
394 
395  //Either Initialize() was not called, or we are in a slave and InitializeLocal() was
396  //not invoked
398  {
399  //create a **thread-local** version of the table. Used only for G4EmCalculator and
400  //Unit Tests
401  fLocalTable = true;
403  }
404 
405  const G4PenelopeCrossSection* theXS =
407  cutEnergy);
408 
409  if (!theXS)
410  {
411  //If we are here, it means that Initialize() was inkoved, but the MaterialTable was
412  //not filled up. This can happen in a UnitTest or via G4EmCalculator
413  if (verboseLevel > 0)
414  {
415  //Issue a G4Exception (warning) only in verbose mode
417  ed << "Unable to retrieve the cross section table for " << theParticle->GetParticleName() <<
418  " in " << material->GetName() << ", cut = " << cutEnergy/keV << " keV " << G4endl;
419  ed << "This can happen only in Unit Tests or via G4EmCalculator" << G4endl;
420  G4Exception("G4PenelopeIonisationModel::ComputeDEDXPerVolume()",
421  "em2038",JustWarning,ed);
422  }
423  //protect file reading via autolock
424  G4AutoLock lock(&PenelopeIonisationModelMutex);
425  theCrossSectionHandler->BuildXSTable(material,cutEnergy,theParticle);
426  lock.unlock();
427  //now it should be ok
428  theXS =
430  material,
431  cutEnergy);
432  }
433 
434  G4double sPowerPerMolecule = 0.0;
435  if (theXS)
436  sPowerPerMolecule = theXS->GetSoftStoppingPower(kineticEnergy);
437 
438 
439  G4double atomDensity = material->GetTotNbOfAtomsPerVolume();
440  G4double atPerMol = oscManager->GetAtomsPerMolecule(material);
441 
442  G4double moleculeDensity = 0.;
443  if (atPerMol)
444  moleculeDensity = atomDensity/atPerMol;
445 
446  G4double sPowerPerVolume = sPowerPerMolecule*moleculeDensity;
447 
448  if (verboseLevel > 2)
449  {
450  G4cout << "G4PenelopeIonisationModel " << G4endl;
451  G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
452  kineticEnergy/keV << " keV = " <<
453  sPowerPerVolume/(keV/mm) << " keV/mm" << G4endl;
454  }
455  return sPowerPerVolume;
456 }
457 
458 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
459 
461  const G4MaterialCutsCouple*)
462 {
464 }
465 
466 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
467 
468 void G4PenelopeIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
469  const G4MaterialCutsCouple* couple,
470  const G4DynamicParticle* aDynamicParticle,
471  G4double cutE, G4double)
472 {
473  // Penelope model v2008 to sample the final state following an hard inelastic interaction.
474  // It makes use of the Generalised Oscillator Strength (GOS) model from
475  // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
476  //
477  // The GOS model is used to calculate the individual cross sections for all
478  // the atomic oscillators coming in the play, taking into account the three
479  // contributions (distant longitudinal collisions, distant transverse collisions and
480  // close collisions). Then the target shell and the interaction channel are
481  // sampled. Final state of the delta-ray (energy, angle) are generated according
482  // to the analytical distributions (dSigma/dW) for the selected interaction
483  // channels.
484  // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because
485  // particles are indistinghusbable), while it is the full initialEnergy for
486  // e+.
487  // The efficiency of the random sampling algorithm (e.g. for close collisions)
488  // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary
489  // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold).
490  // Differential cross sections have a different form for e+ and e-.
491  //
492  // WARNING: The model provides an _average_ description of inelastic collisions.
493  // Anyway, the energy spectrum associated to distant excitations of a given
494  // atomic shell is approximated as a single resonance. The simulated energy spectra
495  // show _unphysical_ narrow peaks at energies that are multiple of the shell
496  // resonance energies. The spurious speaks are automatically smoothed out after
497  // multiple inelastic collisions.
498  //
499  // The model determines also the original shell from which the delta-ray is expelled,
500  // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
501  //
502  // Fermi density correction is calculated analytically according to
503  // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
504 
505  if (verboseLevel > 3)
506  G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl;
507 
508  G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
509  const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
510 
511  if (kineticEnergy0 <= fIntrinsicLowEnergyLimit)
512  {
515  return ;
516  }
517 
518  const G4Material* material = couple->GetMaterial();
520 
521  G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
522 
523  //Initialise final-state variables. The proper values will be set by the methods
524  // SampleFinalStateElectron() and SampleFinalStatePositron()
525  kineticEnergy1=kineticEnergy0;
526  cosThetaPrimary=1.0;
527  energySecondary=0.0;
528  cosThetaSecondary=1.0;
529  targetOscillator = -1;
530 
531  if (theParticle == G4Electron::Electron())
532  SampleFinalStateElectron(material,cutE,kineticEnergy0);
533  else if (theParticle == G4Positron::Positron())
534  SampleFinalStatePositron(material,cutE,kineticEnergy0);
535  else
536  {
538  ed << "Invalid particle " << theParticle->GetParticleName() << G4endl;
539  G4Exception("G4PenelopeIonisationModel::SamplingSecondaries()",
540  "em0001",FatalException,ed);
541 
542  }
543  if (energySecondary == 0) return;
544 
545  if (verboseLevel > 3)
546  {
547  G4cout << "G4PenelopeIonisationModel::SamplingSecondaries() for " <<
548  theParticle->GetParticleName() << G4endl;
549  G4cout << "Final eKin = " << kineticEnergy1 << " keV" << G4endl;
550  G4cout << "Final cosTheta = " << cosThetaPrimary << G4endl;
551  G4cout << "Delta-ray eKin = " << energySecondary << " keV" << G4endl;
552  G4cout << "Delta-ray cosTheta = " << cosThetaSecondary << G4endl;
553  G4cout << "Oscillator: " << targetOscillator << G4endl;
554  }
555 
556  //Update the primary particle
557  G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
558  G4double phiPrimary = twopi * G4UniformRand();
559  G4double dirx = sint * std::cos(phiPrimary);
560  G4double diry = sint * std::sin(phiPrimary);
561  G4double dirz = cosThetaPrimary;
562 
563  G4ThreeVector electronDirection1(dirx,diry,dirz);
564  electronDirection1.rotateUz(particleDirection0);
565 
566  if (kineticEnergy1 > 0)
567  {
568  fParticleChange->ProposeMomentumDirection(electronDirection1);
570  }
571  else
573 
574 
575  //Generate the delta ray
576  G4double ionEnergyInPenelopeDatabase =
577  (*theTable)[targetOscillator]->GetIonisationEnergy();
578 
579  //Now, try to handle fluorescence
580  //Notice: merged levels are indicated with Z=0 and flag=30
581  G4int shFlag = (*theTable)[targetOscillator]->GetShellFlag();
582  G4int Z = (G4int) (*theTable)[targetOscillator]->GetParentZ();
583 
584  //initialize here, then check photons created by Atomic-Deexcitation, and the final state e-
587  //G4int shellId = 0;
588 
589  const G4AtomicShell* shell = 0;
590  //Real level
591  if (Z > 0 && shFlag<30)
592  {
593  shell = transitionManager->Shell(Z,shFlag-1);
594  bindingEnergy = shell->BindingEnergy();
595  //shellId = shell->ShellId();
596  }
597 
598  //correct the energySecondary to account for the fact that the Penelope
599  //database of ionisation energies is in general (slightly) different
600  //from the fluorescence database used in Geant4.
601  energySecondary += ionEnergyInPenelopeDatabase-bindingEnergy;
602 
603  G4double localEnergyDeposit = bindingEnergy;
604  //testing purposes only
605  G4double energyInFluorescence = 0;
606  G4double energyInAuger = 0;
607 
608  if (energySecondary < 0)
609  {
610  //It means that there was some problem/mismatch between the two databases.
611  //In this case, the available energy is ok to excite the level according
612  //to the Penelope database, but not according to the Geant4 database
613  //Full residual energy is deposited locally
614  localEnergyDeposit += energySecondary;
615  energySecondary = 0.0;
616  }
617 
618  //Notice: shell might be nullptr (invalid!) if shFlag=30. Must be protected
619  //Disable the built-in de-excitation of the PIXE flag is active. In this
620  //case, the PIXE interface takes care (statistically) of producing the
621  //de-excitation.
622  //Now, take care of fluorescence, if required
623  if (fAtomDeexcitation && !fPIXEflag && shell)
624  {
625  G4int index = couple->GetIndex();
627  {
628  size_t nBefore = fvect->size();
629  fAtomDeexcitation->GenerateParticles(fvect,shell,Z,index);
630  size_t nAfter = fvect->size();
631 
632  if (nAfter > nBefore) //actual production of fluorescence
633  {
634  for (size_t j=nBefore;j<nAfter;j++) //loop on products
635  {
636  G4double itsEnergy = ((*fvect)[j])->GetKineticEnergy();
637  localEnergyDeposit -= itsEnergy;
638  if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
639  energyInFluorescence += itsEnergy;
640  else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
641  energyInAuger += itsEnergy;
642  }
643  }
644  }
645  }
646 
647  // Generate the delta ray --> to be done only if above cut
648  if (energySecondary > cutE)
649  {
651  G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
652  G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron
653  G4double xEl = sinThetaE * std::cos(phiEl);
654  G4double yEl = sinThetaE * std::sin(phiEl);
656  G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
657  eDirection.rotateUz(particleDirection0);
658  electron = new G4DynamicParticle (G4Electron::Electron(),
659  eDirection,energySecondary) ;
660  fvect->push_back(electron);
661  }
662  else
663  {
664  localEnergyDeposit += energySecondary;
665  energySecondary = 0;
666  }
667 
668  if (localEnergyDeposit < 0)
669  {
670  G4cout << "WARNING-"
671  << "G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
672  << G4endl;
673  localEnergyDeposit=0.;
674  }
675  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
676 
677  if (verboseLevel > 1)
678  {
679  G4cout << "-----------------------------------------------------------" << G4endl;
680  G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
681  G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
682  G4cout << "-----------------------------------------------------------" << G4endl;
683  G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
684  G4cout << "Delta ray " << energySecondary/keV << " keV" << G4endl;
685  if (energyInFluorescence)
686  G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
687  if (energyInAuger)
688  G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
689  G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
690  G4cout << "Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+
691  localEnergyDeposit+energyInAuger)/keV <<
692  " keV" << G4endl;
693  G4cout << "-----------------------------------------------------------" << G4endl;
694  }
695 
696  if (verboseLevel > 0)
697  {
698  G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
699  localEnergyDeposit+energyInAuger-kineticEnergy0);
700  if (energyDiff > 0.05*keV)
701  G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<
702  (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit+energyInAuger)/keV <<
703  " keV (final) vs. " <<
704  kineticEnergy0/keV << " keV (initial)" << G4endl;
705  }
706 
707 }
708 
709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
711  G4double cutEnergy,
712  G4double kineticEnergy)
713 {
714  // This method sets the final ionisation parameters
715  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
716  // energySecondary, cosThetaSecondary (= info of delta-ray)
717  // targetOscillator (= ionised oscillator)
718  //
719  // The method implements SUBROUTINE EINa of Penelope
720  //
721 
723  size_t numberOfOscillators = theTable->size();
724  const G4PenelopeCrossSection* theXS =
726  cutEnergy);
727  G4double delta = theCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
728 
729  // Selection of the active oscillator
730  G4double TST = G4UniformRand();
731  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
732  G4double XSsum = 0.;
733 
734  for (size_t i=0;i<numberOfOscillators-1;i++)
735  {
736  /* testing purposes
737  G4cout << "sampling: " << i << " " << XSsum << " " << TST << " " <<
738  theXS->GetShellCrossSection(i,kineticEnergy) << " " <<
739  theXS->GetNormalizedShellCrossSection(i,kineticEnergy) << " " <<
740  mat->GetName() << G4endl;
741  */
742  XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
743 
744  if (XSsum > TST)
745  {
746  targetOscillator = (G4int) i;
747  break;
748  }
749  }
750 
751 
752  if (verboseLevel > 3)
753  {
754  G4cout << "SampleFinalStateElectron: sampled oscillator #" << targetOscillator << "." << G4endl;
755  G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV <<
756  " eV " << G4endl;
757  G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV << " eV "
758  << G4endl;
759  }
760  //Constants
761  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
762  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
763  G4double gam2 = gam*gam;
764  G4double beta2 = (gam2-1.0)/gam2;
765  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
766 
767  //Partial cross section of the active oscillator
768  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
769  G4double invResEne = 1.0/resEne;
770  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
771  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
772  G4double XHDL = 0.;
773  G4double XHDT = 0.;
774  G4double QM = 0.;
775  G4double cps = 0.;
776  G4double cp = 0.;
777 
778  //Distant excitations
779  if (resEne > cutEnergy && resEne < kineticEnergy)
780  {
781  cps = kineticEnergy*rb;
782  cp = std::sqrt(cps);
783  G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
784  if (resEne > 1.0e-6*kineticEnergy)
785  {
786  G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
787  QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
788  }
789  else
790  {
791  QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
792  QM *= (1.0-QM*0.5/electron_mass_c2);
793  }
794  if (QM < cutoffEne)
795  {
796  XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
797  *invResEne;
798  XHDT = XHDT0*invResEne;
799  }
800  else
801  {
802  QM = cutoffEne;
803  XHDL = 0.;
804  XHDT = 0.;
805  }
806  }
807  else
808  {
809  QM = cutoffEne;
810  cps = 0.;
811  cp = 0.;
812  XHDL = 0.;
813  XHDT = 0.;
814  }
815 
816  //Close collisions
817  G4double EE = kineticEnergy + ionEne;
818  G4double wmaxc = 0.5*EE;
819  G4double wcl = std::max(cutEnergy,cutoffEne);
820  G4double rcl = wcl/EE;
821  G4double XHC = 0.;
822  if (wcl < wmaxc)
823  {
824  G4double rl1 = 1.0-rcl;
825  G4double rrl1 = 1.0/rl1;
826  XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
827  (1.0-amol)*std::log(rcl*rrl1))/EE;
828  }
829 
830  //Total cross section per molecule for the active shell, in cm2
831  G4double XHTOT = XHC + XHDL + XHDT;
832 
833  //very small cross section, do nothing
834  if (XHTOT < 1.e-14*barn)
835  {
836  kineticEnergy1=kineticEnergy;
837  cosThetaPrimary=1.0;
838  energySecondary=0.0;
839  cosThetaSecondary=1.0;
840  targetOscillator = numberOfOscillators-1;
841  return;
842  }
843 
844  //decide which kind of interaction we'll have
845  TST = XHTOT*G4UniformRand();
846 
847 
848  // Hard close collision
849  G4double TS1 = XHC;
850 
851  if (TST < TS1)
852  {
853  G4double A = 5.0*amol;
854  G4double ARCL = A*0.5*rcl;
855  G4double rk=0.;
856  G4bool loopAgain = false;
857  do
858  {
859  loopAgain = false;
860  G4double fb = (1.0+ARCL)*G4UniformRand();
861  if (fb < 1)
862  rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
863  else
864  rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
865  G4double rk2 = rk*rk;
866  G4double rkf = rk/(1.0-rk);
867  G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
868  if (G4UniformRand()*(1.0+A*rk2) > phi)
869  loopAgain = true;
870  }while(loopAgain);
871  //energy and scattering angle (primary electron)
872  G4double deltaE = rk*EE;
873 
874  kineticEnergy1 = kineticEnergy - deltaE;
875  cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
876  //energy and scattering angle of the delta ray
877  energySecondary = deltaE - ionEne; //subtract ionisation energy
878  cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
879  if (verboseLevel > 3)
880  G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl;
881  return;
882  }
883 
884  //Hard distant longitudinal collisions
885  TS1 += XHDL;
886  G4double deltaE = resEne;
887  kineticEnergy1 = kineticEnergy - deltaE;
888 
889  if (TST < TS1)
890  {
891  G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
892  G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
893  - (QS*0.5/electron_mass_c2));
894  G4double QTREV = Q*(Q+2.0*electron_mass_c2);
895  G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
896  cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
897  if (cosThetaPrimary > 1.)
898  cosThetaPrimary = 1.0;
899  //energy and emission angle of the delta ray
900  energySecondary = deltaE - ionEne;
901  cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
902  if (cosThetaSecondary > 1.0)
903  cosThetaSecondary = 1.0;
904  if (verboseLevel > 3)
905  G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl;
906  return;
907  }
908 
909  //Hard distant transverse collisions
910  cosThetaPrimary = 1.0;
911  //energy and emission angle of the delta ray
912  energySecondary = deltaE - ionEne;
913  cosThetaSecondary = 0.5;
914  if (verboseLevel > 3)
915  G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl;
916 
917  return;
918 }
919 
920 
921 
922 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
924  G4double cutEnergy,
925  G4double kineticEnergy)
926 {
927  // This method sets the final ionisation parameters
928  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
929  // energySecondary, cosThetaSecondary (= info of delta-ray)
930  // targetOscillator (= ionised oscillator)
931  //
932  // The method implements SUBROUTINE PINa of Penelope
933  //
934 
936  size_t numberOfOscillators = theTable->size();
937  const G4PenelopeCrossSection* theXS =
939  cutEnergy);
940  G4double delta = theCrossSectionHandler->GetDensityCorrection(mat,kineticEnergy);
941 
942  // Selection of the active oscillator
943  G4double TST = G4UniformRand();
944  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
945  G4double XSsum = 0.;
946  for (size_t i=0;i<numberOfOscillators-1;i++)
947  {
948  XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
949  if (XSsum > TST)
950  {
951  targetOscillator = (G4int) i;
952  break;
953  }
954  }
955 
956  if (verboseLevel > 3)
957  {
958  G4cout << "SampleFinalStatePositron: sampled oscillator #" << targetOscillator << "." << G4endl;
959  G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV
960  << " eV " << G4endl;
961  G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV
962  << " eV " << G4endl;
963  }
964 
965  //Constants
966  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
967  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
968  G4double gam2 = gam*gam;
969  G4double beta2 = (gam2-1.0)/gam2;
970  G4double g12 = (gam+1.0)*(gam+1.0);
971  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
972  //Bhabha coefficients
973  G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
974  G4double bha2 = amol*(3.0+1.0/g12);
975  G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
976  G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
977 
978  //
979  //Partial cross section of the active oscillator
980  //
981  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
982  G4double invResEne = 1.0/resEne;
983  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
984  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
985 
986  G4double XHDL = 0.;
987  G4double XHDT = 0.;
988  G4double QM = 0.;
989  G4double cps = 0.;
990  G4double cp = 0.;
991 
992  //Distant excitations XS (same as for electrons)
993  if (resEne > cutEnergy && resEne < kineticEnergy)
994  {
995  cps = kineticEnergy*rb;
996  cp = std::sqrt(cps);
997  G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
998  if (resEne > 1.0e-6*kineticEnergy)
999  {
1000  G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
1001  QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
1002  }
1003  else
1004  {
1005  QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
1006  QM *= (1.0-QM*0.5/electron_mass_c2);
1007  }
1008  if (QM < cutoffEne)
1009  {
1010  XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
1011  *invResEne;
1012  XHDT = XHDT0*invResEne;
1013  }
1014  else
1015  {
1016  QM = cutoffEne;
1017  XHDL = 0.;
1018  XHDT = 0.;
1019  }
1020  }
1021  else
1022  {
1023  QM = cutoffEne;
1024  cps = 0.;
1025  cp = 0.;
1026  XHDL = 0.;
1027  XHDT = 0.;
1028  }
1029  //Close collisions (Bhabha)
1030  G4double wmaxc = kineticEnergy;
1031  G4double wcl = std::max(cutEnergy,cutoffEne);
1032  G4double rcl = wcl/kineticEnergy;
1033  G4double XHC = 0.;
1034  if (wcl < wmaxc)
1035  {
1036  G4double rl1 = 1.0-rcl;
1037  XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
1038  + (bha3/2.0)*(rcl*rcl-1.0)
1039  + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1040  }
1041 
1042  //Total cross section per molecule for the active shell, in cm2
1043  G4double XHTOT = XHC + XHDL + XHDT;
1044 
1045  //very small cross section, do nothing
1046  if (XHTOT < 1.e-14*barn)
1047  {
1048  kineticEnergy1=kineticEnergy;
1049  cosThetaPrimary=1.0;
1050  energySecondary=0.0;
1051  cosThetaSecondary=1.0;
1052  targetOscillator = numberOfOscillators-1;
1053  return;
1054  }
1055 
1056  //decide which kind of interaction we'll have
1057  TST = XHTOT*G4UniformRand();
1058 
1059  // Hard close collision
1060  G4double TS1 = XHC;
1061  if (TST < TS1)
1062  {
1063  G4double rl1 = 1.0-rcl;
1064  G4double rk=0.;
1065  G4bool loopAgain = false;
1066  do
1067  {
1068  loopAgain = false;
1069  rk = rcl/(1.0-G4UniformRand()*rl1);
1070  G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1071  if (G4UniformRand() > phi)
1072  loopAgain = true;
1073  }while(loopAgain);
1074  //energy and scattering angle (primary electron)
1075  G4double deltaE = rk*kineticEnergy;
1076  kineticEnergy1 = kineticEnergy - deltaE;
1077  cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1078  //energy and scattering angle of the delta ray
1079  energySecondary = deltaE - ionEne; //subtract ionisation energy
1080  cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
1081  if (verboseLevel > 3)
1082  G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl;
1083  return;
1084  }
1085 
1086  //Hard distant longitudinal collisions
1087  TS1 += XHDL;
1088  G4double deltaE = resEne;
1089  kineticEnergy1 = kineticEnergy - deltaE;
1090  if (TST < TS1)
1091  {
1092  G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
1093  G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
1094  - (QS*0.5/electron_mass_c2));
1095  G4double QTREV = Q*(Q+2.0*electron_mass_c2);
1096  G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
1097  cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1098  if (cosThetaPrimary > 1.)
1099  cosThetaPrimary = 1.0;
1100  //energy and emission angle of the delta ray
1101  energySecondary = deltaE - ionEne;
1102  cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1103  if (cosThetaSecondary > 1.0)
1104  cosThetaSecondary = 1.0;
1105  if (verboseLevel > 3)
1106  G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl;
1107  return;
1108  }
1109 
1110  //Hard distant transverse collisions
1111  cosThetaPrimary = 1.0;
1112  //energy and emission angle of the delta ray
1113  energySecondary = deltaE - ionEne;
1114  cosThetaSecondary = 0.5;
1115 
1116  if (verboseLevel > 3)
1117  G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl;
1118 
1119  return;
1120 }
1121 
1122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1123 
1125 {
1126  if(!fParticle) {
1127  fParticle = p;
1128  }
1129 }
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4VAtomDeexcitation * fAtomDeexcitation
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
static constexpr double mm
Definition: G4SIunits.hh:115
G4double GetSoftStoppingPower(G4double energy) const
Returns the total stopping power due to soft collisions.
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
void SetVerboseLevel(G4int vl)
Setter for the verbosity level.
CLHEP::Hep3Vector G4ThreeVector
void SampleFinalStateElectron(const G4Material *, G4double cutEnergy, G4double kineticEnergy)
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
virtual G4double CrossSectionPerVolume(const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
G4bool IsPIXEActive() const
G4PenelopeOscillatorManager * oscManager
const G4String & GetName() const
Definition: G4Material.hh:178
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:411
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
static double Q[]
G4bool IsMaster() const
Definition: G4VEmModel.hh:717
G4ParticleDefinition * GetDefinition() const
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
Returns the table of cross sections for the given particle, given material and given cut as a G4Penel...
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *)
G4double BindingEnergy() const
int G4int
Definition: G4Types.hh:78
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:724
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4PenelopeIonisationXSHandler * theCrossSectionHandler
G4double GetTotalCrossSection(G4double energy) const
Returns total cross section at the given energy.
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double eV
Definition: G4SIunits.hh:215
G4ParticleChangeForLoss * fParticleChange
static G4PenelopeOscillatorManager * GetOscillatorManager()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void SetParticle(const G4ParticleDefinition *)
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
static G4ProductionCutsTable * GetProductionCutsTable()
virtual G4double MinEnergyCut(const G4ParticleDefinition *, const G4MaterialCutsCouple *)
G4int G4Mutex
Definition: G4Threading.hh:173
static G4Positron * Positron()
Definition: G4Positron.cc:94
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
static G4EmParameters * Instance()
G4double GetDensityCorrection(const G4Material *, const G4double energy) const
Returns the density coeection for the material at the given energy.
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
G4VAtomDeexcitation * AtomDeexcitation()
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
void BuildXSTable(const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
This can be inkoved only by the master.
double G4double
Definition: G4Types.hh:76
G4double GetNormalizedShellCrossSection(size_t shellID, G4double energy) const
Returns the hard cross section for the given shell (normalized to 1)
const G4ParticleDefinition * fParticle
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:780
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
static constexpr double barn
Definition: G4SIunits.hh:105
G4double GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
G4ThreeVector G4ParticleMomentum
G4double bindingEnergy(G4int A, G4int Z)
G4PenelopeIonisationModel(const G4ParticleDefinition *p=0, const G4String &processName="PenIoni")
void SampleFinalStatePositron(const G4Material *, G4double cutEnergy, G4double kineticEnergy)
static G4AtomicTransitionManager * Instance()
static constexpr double keV
Definition: G4SIunits.hh:216
const G4String & PIXEElectronCrossSectionModel()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
static const G4double cp
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
const G4Material * GetMaterial() const