Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeIonisationModel Class Reference

#include <G4PenelopeIonisationModel.hh>

Inheritance diagram for G4PenelopeIonisationModel:
Collaboration diagram for G4PenelopeIonisationModel:

Public Member Functions

 G4PenelopeIonisationModel (const G4ParticleDefinition *p=0, const G4String &processName="PenIoni")
 
virtual ~G4PenelopeIonisationModel ()
 
virtual void Initialise (const G4ParticleDefinition *, const G4DataVector &)
 
virtual void InitialiseLocal (const G4ParticleDefinition *, G4VEmModel *)
 
virtual G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double, G4double, G4double, G4double, G4double)
 
virtual G4double CrossSectionPerVolume (const G4Material *material, const G4ParticleDefinition *theParticle, G4double kineticEnergy, G4double cutEnergy, G4double maxEnergy=DBL_MAX)
 
virtual void SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
 
virtual G4double ComputeDEDXPerVolume (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
 
virtual G4double MinEnergyCut (const G4ParticleDefinition *, const G4MaterialCutsCouple *)
 
void SetVerbosityLevel (G4int lev)
 
G4int GetVerbosityLevel ()
 
- Public Member Functions inherited from G4VEmModel
 G4VEmModel (const G4String &nam)
 
virtual ~G4VEmModel ()
 
virtual void InitialiseForMaterial (const G4ParticleDefinition *, const G4Material *)
 
virtual void InitialiseForElement (const G4ParticleDefinition *, G4int Z)
 
virtual G4double GetPartialCrossSection (const G4Material *, G4int level, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double ComputeCrossSectionPerShell (const G4ParticleDefinition *, G4int Z, G4int shellIdx, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
virtual G4double ChargeSquareRatio (const G4Track &)
 
virtual G4double GetChargeSquareRatio (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual G4double GetParticleCharge (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void StartTracking (G4Track *)
 
virtual void CorrectionsAlongStep (const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double &eloss, G4double &niel, G4double length)
 
virtual G4double Value (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy)
 
virtual G4double MinPrimaryEnergy (const G4Material *, const G4ParticleDefinition *, G4double cut=0.0)
 
virtual void SetupForMaterial (const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
 
virtual void DefineForRegion (const G4Region *)
 
virtual void ModelDescription (std::ostream &outFile) const
 
void InitialiseElementSelectors (const G4ParticleDefinition *, const G4DataVector &)
 
std::vector
< G4EmElementSelector * > * 
GetElementSelectors ()
 
void SetElementSelectors (std::vector< G4EmElementSelector * > *)
 
virtual G4double ComputeDEDX (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=DBL_MAX)
 
G4double CrossSection (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeMeanFreePath (const G4ParticleDefinition *, G4double kineticEnergy, const G4Material *, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4double ComputeCrossSectionPerAtom (const G4ParticleDefinition *, const G4Element *, G4double kinEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectIsotopeNumber (const G4Element *)
 
const G4ElementSelectRandomAtom (const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
const G4ElementSelectRandomAtom (const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
 
G4int SelectRandomAtomNumber (const G4Material *)
 
void SetParticleChange (G4VParticleChange *, G4VEmFluctuationModel *f=nullptr)
 
void SetCrossSectionTable (G4PhysicsTable *, G4bool isLocal)
 
G4ElementDataGetElementData ()
 
G4PhysicsTableGetCrossSectionTable ()
 
G4VEmFluctuationModelGetModelOfFluctuations ()
 
G4VEmAngularDistributionGetAngularDistribution ()
 
void SetAngularDistribution (G4VEmAngularDistribution *)
 
G4double HighEnergyLimit () const
 
G4double LowEnergyLimit () const
 
G4double HighEnergyActivationLimit () const
 
G4double LowEnergyActivationLimit () const
 
G4double PolarAngleLimit () const
 
G4double SecondaryThreshold () const
 
G4bool LPMFlag () const
 
G4bool DeexcitationFlag () const
 
G4bool ForceBuildTableFlag () const
 
G4bool UseAngularGeneratorFlag () const
 
void SetAngularGeneratorFlag (G4bool)
 
void SetHighEnergyLimit (G4double)
 
void SetLowEnergyLimit (G4double)
 
void SetActivationHighEnergyLimit (G4double)
 
void SetActivationLowEnergyLimit (G4double)
 
G4bool IsActive (G4double kinEnergy)
 
void SetPolarAngleLimit (G4double)
 
void SetSecondaryThreshold (G4double)
 
void SetLPMFlag (G4bool val)
 
void SetDeexcitationFlag (G4bool val)
 
void SetForceBuildTable (G4bool val)
 
void SetFluctuationFlag (G4bool val)
 
void SetMasterThread (G4bool val)
 
G4bool IsMaster () const
 
G4double MaxSecondaryKinEnergy (const G4DynamicParticle *dynParticle)
 
const G4StringGetName () const
 
void SetCurrentCouple (const G4MaterialCutsCouple *)
 
const G4ElementGetCurrentElement () const
 
const G4IsotopeGetCurrentIsotope () const
 
G4bool IsLocked () const
 
void SetLocked (G4bool)
 

Protected Attributes

G4ParticleChangeForLossfParticleChange
 
const G4ParticleDefinitionfParticle
 
- Protected Attributes inherited from G4VEmModel
G4ElementDatafElementData
 
G4VParticleChangepParticleChange
 
G4PhysicsTablexSectionTable
 
const std::vector< G4double > * theDensityFactor
 
const std::vector< G4int > * theDensityIdx
 
size_t idxTable
 
G4bool lossFlucFlag
 

Additional Inherited Members

- Protected Member Functions inherited from G4VEmModel
G4ParticleChangeForLossGetParticleChangeForLoss ()
 
G4ParticleChangeForGammaGetParticleChangeForGamma ()
 
virtual G4double MaxSecondaryEnergy (const G4ParticleDefinition *, G4double kineticEnergy)
 
const G4MaterialCutsCoupleCurrentCouple () const
 
void SetCurrentElement (const G4Element *)
 
- Static Protected Attributes inherited from G4VEmModel
static const G4double inveplus = 1.0/CLHEP::eplus
 

Detailed Description

Definition at line 66 of file G4PenelopeIonisationModel.hh.

Constructor & Destructor Documentation

G4PenelopeIonisationModel::G4PenelopeIonisationModel ( const G4ParticleDefinition p = 0,
const G4String processName = "PenIoni" 
)

Definition at line 73 of file G4PenelopeIonisationModel.cc.

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 {
81  fIntrinsicLowEnergyLimit = 100.0*eV;
82  fIntrinsicHighEnergyLimit = 100.0*GeV;
83  // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
84  SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
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 }
G4VEmModel(const G4String &nam)
Definition: G4VEmModel.cc:68
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
static constexpr double eV
Definition: G4SIunits.hh:215
G4ParticleChangeForLoss * fParticleChange
static G4PenelopeOscillatorManager * GetOscillatorManager()
static constexpr double GeV
Definition: G4SIunits.hh:217
const G4ParticleDefinition * fParticle
void SetDeexcitationFlag(G4bool val)
Definition: G4VEmModel.hh:788

Here is the call graph for this function:

G4PenelopeIonisationModel::~G4PenelopeIonisationModel ( )
virtual

Definition at line 108 of file G4PenelopeIonisationModel.cc.

109 {
110  if (IsMaster() || fLocalTable)
111  {
112  if (theCrossSectionHandler)
113  delete theCrossSectionHandler;
114  }
115 }
G4bool IsMaster() const
Definition: G4VEmModel.hh:725

Here is the call graph for this function:

Member Function Documentation

G4double G4PenelopeIonisationModel::ComputeCrossSectionPerAtom ( const G4ParticleDefinition ,
G4double  ,
G4double  ,
G4double  ,
G4double  ,
G4double   
)
virtual

Reimplemented from G4VEmModel.

Definition at line 356 of file G4PenelopeIonisationModel.cc.

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 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume ( const G4Material material,
const G4ParticleDefinition theParticle,
G4double  kineticEnergy,
G4double  cutEnergy 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 372 of file G4PenelopeIonisationModel.cc.

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
397  if (!theCrossSectionHandler)
398  {
399  //create a **thread-local** version of the table. Used only for G4EmCalculator and
400  //Unit Tests
401  fLocalTable = true;
402  theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);
403  }
404 
405  const G4PenelopeCrossSection* theXS =
406  theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,material,
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 =
429  theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
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 }
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
const G4String & GetName() const
Definition: G4Material.hh:178
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
#define G4endl
Definition: G4ios.hh:61
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 GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4double G4PenelopeIonisationModel::CrossSectionPerVolume ( const G4Material material,
const G4ParticleDefinition theParticle,
G4double  kineticEnergy,
G4double  cutEnergy,
G4double  maxEnergy = DBL_MAX 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 251 of file G4PenelopeIonisationModel.cc.

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
284  if (!theCrossSectionHandler)
285  {
286  //create a **thread-local** version of the table. Used only for G4EmCalculator and
287  //Unit Tests
288  fLocalTable = true;
289  theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);
290  }
291 
292  const G4PenelopeCrossSection* theXS =
293  theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
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 =
316  theCrossSectionHandler->GetCrossSectionTableForCouple(theParticle,
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 }
static constexpr double mm
Definition: G4SIunits.hh:115
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
const G4String & GetName() const
Definition: G4Material.hh:178
virtual void SetupForMaterial(const G4ParticleDefinition *, const G4Material *, G4double kineticEnergy)
Definition: G4VEmModel.cc:411
G4double GetHardCrossSection(G4double energy) const
Returns hard cross section at the given energy.
const G4PenelopeCrossSection * GetCrossSectionTableForCouple(const G4ParticleDefinition *, const G4Material *, const G4double cut) const
const G4String & GetParticleName() const
G4GLOB_DLL std::ostream G4cout
G4double GetTotalCrossSection(G4double energy) const
Returns total cross section at the given energy.
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:209
G4double energy(const ThreeVector &p, const G4double m)
#define G4endl
Definition: G4ios.hh:61
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 GetAtomsPerMolecule(const G4Material *)
Returns the total number of atoms per molecule.
static constexpr double keV
Definition: G4SIunits.hh:216

Here is the call graph for this function:

G4int G4PenelopeIonisationModel::GetVerbosityLevel ( )
inline

Definition at line 113 of file G4PenelopeIonisationModel.hh.

113 {return verboseLevel;};
void G4PenelopeIonisationModel::Initialise ( const G4ParticleDefinition particle,
const G4DataVector theCuts 
)
virtual

Implements G4VEmModel.

Definition at line 119 of file G4PenelopeIonisationModel.cc.

121 {
122  if (verboseLevel > 3)
123  G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
124 
125  fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
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)
137  fPIXEflag = fAtomDeexcitation->IsPIXEActive();
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
178  if (theCrossSectionHandler)
179  {
180  delete theCrossSectionHandler;
181  theCrossSectionHandler = 0;
182  }
183  theCrossSectionHandler = new G4PenelopeIonisationXSHandler(nBins);
184  theCrossSectionHandler->SetVerboseLevel(verboseLevel);
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 }
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:643
static G4LossTableManager * Instance()
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
void SetVerboseLevel(G4int vl)
Setter for the verbosity level.
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:636
G4bool IsPIXEActive() const
G4bool IsMaster() const
Definition: G4VEmModel.hh:725
G4GLOB_DLL std::ostream G4cout
G4ParticleChangeForLoss * fParticleChange
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static G4EmParameters * Instance()
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
void BuildXSTable(const G4Material *, G4double cut, const G4ParticleDefinition *, G4bool isMaster=true)
This can be inkoved only by the master.
const G4ParticleDefinition * fParticle
static constexpr double keV
Definition: G4SIunits.hh:216
const G4String & PIXEElectronCrossSectionModel()
const G4Material * GetMaterial() const

Here is the call graph for this function:

void G4PenelopeIonisationModel::InitialiseLocal ( const G4ParticleDefinition part,
G4VEmModel masterModel 
)
virtual

Reimplemented from G4VEmModel.

Definition at line 219 of file G4PenelopeIonisationModel.cc.

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
236  theCrossSectionHandler = theModel->theCrossSectionHandler;
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 }
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
const G4ParticleDefinition * fParticle
G4double G4PenelopeIonisationModel::MinEnergyCut ( const G4ParticleDefinition ,
const G4MaterialCutsCouple  
)
virtual

Reimplemented from G4VEmModel.

Definition at line 460 of file G4PenelopeIonisationModel.cc.

462 {
463  return fIntrinsicLowEnergyLimit;
464 }
void G4PenelopeIonisationModel::SampleSecondaries ( std::vector< G4DynamicParticle * > *  fvect,
const G4MaterialCutsCouple couple,
const G4DynamicParticle aDynamicParticle,
G4double  tmin,
G4double  maxEnergy 
)
virtual

Implements G4VEmModel.

Definition at line 468 of file G4PenelopeIonisationModel.cc.

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();
519  const G4PenelopeOscillatorTable* theTable = oscManager->GetOscillatorTableIonisation(material);
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);
569  fParticleChange->SetProposedKineticEnergy(kineticEnergy1);
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();
626  if (fAtomDeexcitation->CheckDeexcitationActiveRegion(index))
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);
655  G4double zEl = cosThetaSecondary;
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 }
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4PenelopeOscillatorTable * GetOscillatorTableIonisation(const G4Material *)
G4bool CheckDeexcitationActiveRegion(G4int coupleIndex)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4ParticleDefinition * GetDefinition() const
G4double BindingEnergy() const
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double twopi
Definition: G4SIunits.hh:76
string material
Definition: eplot.py:19
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double eV
Definition: G4SIunits.hh:215
G4ParticleChangeForLoss * fParticleChange
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4Positron * Positron()
Definition: G4Positron.cc:94
std::vector< G4PenelopeOscillator * > G4PenelopeOscillatorTable
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static constexpr double pi
Definition: G4SIunits.hh:75
double G4double
Definition: G4Types.hh:76
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
G4double bindingEnergy(G4int A, G4int Z)
static G4AtomicTransitionManager * Instance()
static constexpr double keV
Definition: G4SIunits.hh:216
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
const G4Material * GetMaterial() const

Here is the call graph for this function:

void G4PenelopeIonisationModel::SetVerbosityLevel ( G4int  lev)
inline

Definition at line 112 of file G4PenelopeIonisationModel.hh.

112 {verboseLevel = lev;};

Member Data Documentation

const G4ParticleDefinition* G4PenelopeIonisationModel::fParticle
protected

Definition at line 117 of file G4PenelopeIonisationModel.hh.

G4ParticleChangeForLoss* G4PenelopeIonisationModel::fParticleChange
protected

Definition at line 113 of file G4PenelopeIonisationModel.hh.


The documentation for this class was generated from the following files: