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