ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4PenelopeIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4PenelopeIonisationModel.cc
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 //
27 // Author: Luciano Pandola
28 //
29 // History:
30 // --------
31 // 27 Jul 2010 L Pandola First complete implementation
32 // 18 Jan 2011 L.Pandola Stricter check on production of sub-treshold delta-rays.
33 // Should never happen now
34 // 01 Feb 2011 L Pandola Suppress fake energy-violation warning when Auger is active.
35 // Make sure that fluorescence/Auger is generated only if
36 // above threshold
37 // 25 May 2011 L Pandola Renamed (make v2008 as default Penelope)
38 // 26 Jan 2012 L Pandola Migration of AtomicDeexcitation to the new interface
39 // 09 Mar 2012 L Pandola Moved the management and calculation of
40 // cross sections to a separate class. Use a different method to
41 // get normalized shell cross sections
42 // 07 Oct 2013 L. Pandola Migration to MT
43 // 23 Jun 2015 L. Pandola Keep track of the PIXE flag, to avoid double-production of
44 // atomic de-excitation (bug #1761)
45 // 29 Aug 2018 L. Pandola Fix bug causing energy non-conservation
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  if (itsEnergy < localEnergyDeposit) // valid secondary, generate it
638  {
639  localEnergyDeposit -= itsEnergy;
640  if (((*fvect)[j])->GetParticleDefinition() == G4Gamma::Definition())
641  energyInFluorescence += itsEnergy;
642  else if (((*fvect)[j])->GetParticleDefinition() == G4Electron::Definition())
643  energyInAuger += itsEnergy;
644  }
645  else //invalid secondary: takes more than the available energy: delete it
646  {
647  delete (*fvect)[j];
648  (*fvect)[j] = nullptr;
649  }
650 
651  }
652  }
653  }
654  }
655 
656  // Generate the delta ray --> to be done only if above cut
657  if (energySecondary > cutE)
658  {
660  G4double sinThetaE = std::sqrt(1.-cosThetaSecondary*cosThetaSecondary);
661  G4double phiEl = phiPrimary+pi; //pi with respect to the primary electron/positron
662  G4double xEl = sinThetaE * std::cos(phiEl);
663  G4double yEl = sinThetaE * std::sin(phiEl);
665  G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
666  eDirection.rotateUz(particleDirection0);
667  electron = new G4DynamicParticle (G4Electron::Electron(),
668  eDirection,energySecondary) ;
669  fvect->push_back(electron);
670  }
671  else
672  {
673  localEnergyDeposit += energySecondary;
674  energySecondary = 0;
675  }
676 
677  if (localEnergyDeposit < 0) //Should not be: issue a G4Exception (warning)
678  {
679  G4Exception("G4PenelopeIonisationModel::SampleSecondaries()",
680  "em2099",JustWarning,"WARNING: Negative local energy deposit");
681  localEnergyDeposit=0.;
682  }
683  fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
684 
685  if (verboseLevel > 1)
686  {
687  G4cout << "-----------------------------------------------------------" << G4endl;
688  G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
689  G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
690  G4cout << "-----------------------------------------------------------" << G4endl;
691  G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
692  G4cout << "Delta ray " << energySecondary/keV << " keV" << G4endl;
693  if (energyInFluorescence)
694  G4cout << "Fluorescence x-rays: " << energyInFluorescence/keV << " keV" << G4endl;
695  if (energyInAuger)
696  G4cout << "Auger electrons: " << energyInAuger/keV << " keV" << G4endl;
697  G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
698  G4cout << "Total final state: " << (energySecondary+energyInFluorescence+kineticEnergy1+
699  localEnergyDeposit+energyInAuger)/keV <<
700  " keV" << G4endl;
701  G4cout << "-----------------------------------------------------------" << G4endl;
702  }
703 
704  if (verboseLevel > 0)
705  {
706  G4double energyDiff = std::fabs(energySecondary+energyInFluorescence+kineticEnergy1+
707  localEnergyDeposit+energyInAuger-kineticEnergy0);
708  if (energyDiff > 0.05*keV)
709  G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<
710  (energySecondary+energyInFluorescence+kineticEnergy1+localEnergyDeposit+energyInAuger)/keV <<
711  " keV (final) vs. " <<
712  kineticEnergy0/keV << " keV (initial)" << G4endl;
713  }
714 
715 }
716 
717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
719  G4double cutEnergy,
720  G4double kineticEnergy)
721 {
722  // This method sets the final ionisation parameters
723  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
724  // energySecondary, cosThetaSecondary (= info of delta-ray)
725  // targetOscillator (= ionised oscillator)
726  //
727  // The method implements SUBROUTINE EINa of Penelope
728  //
729 
731  size_t numberOfOscillators = theTable->size();
732  const G4PenelopeCrossSection* theXS =
734  cutEnergy);
736 
737  // Selection of the active oscillator
738  G4double TST = G4UniformRand();
739  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
740  G4double XSsum = 0.;
741 
742  for (size_t i=0;i<numberOfOscillators-1;i++)
743  {
744  /* testing purposes
745  G4cout << "sampling: " << i << " " << XSsum << " " << TST << " " <<
746  theXS->GetShellCrossSection(i,kineticEnergy) << " " <<
747  theXS->GetNormalizedShellCrossSection(i,kineticEnergy) << " " <<
748  mat->GetName() << G4endl;
749  */
750  XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
751 
752  if (XSsum > TST)
753  {
754  targetOscillator = (G4int) i;
755  break;
756  }
757  }
758 
759 
760  if (verboseLevel > 3)
761  {
762  G4cout << "SampleFinalStateElectron: sampled oscillator #" << targetOscillator << "." << G4endl;
763  G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV <<
764  " eV " << G4endl;
765  G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV << " eV "
766  << G4endl;
767  }
768  //Constants
769  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
770  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
771  G4double gam2 = gam*gam;
772  G4double beta2 = (gam2-1.0)/gam2;
773  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
774 
775  //Partial cross section of the active oscillator
776  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
777  G4double invResEne = 1.0/resEne;
778  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
779  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
780  G4double XHDL = 0.;
781  G4double XHDT = 0.;
782  G4double QM = 0.;
783  G4double cps = 0.;
784  G4double cp = 0.;
785 
786  //Distant excitations
787  if (resEne > cutEnergy && resEne < kineticEnergy)
788  {
789  cps = kineticEnergy*rb;
790  cp = std::sqrt(cps);
791  G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
792  if (resEne > 1.0e-6*kineticEnergy)
793  {
794  G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
795  QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
796  }
797  else
798  {
799  QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
800  QM *= (1.0-QM*0.5/electron_mass_c2);
801  }
802  if (QM < cutoffEne)
803  {
804  XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
805  *invResEne;
806  XHDT = XHDT0*invResEne;
807  }
808  else
809  {
810  QM = cutoffEne;
811  XHDL = 0.;
812  XHDT = 0.;
813  }
814  }
815  else
816  {
817  QM = cutoffEne;
818  cps = 0.;
819  cp = 0.;
820  XHDL = 0.;
821  XHDT = 0.;
822  }
823 
824  //Close collisions
825  G4double EE = kineticEnergy + ionEne;
826  G4double wmaxc = 0.5*EE;
827  G4double wcl = std::max(cutEnergy,cutoffEne);
828  G4double rcl = wcl/EE;
829  G4double XHC = 0.;
830  if (wcl < wmaxc)
831  {
832  G4double rl1 = 1.0-rcl;
833  G4double rrl1 = 1.0/rl1;
834  XHC = (amol*(0.5-rcl)+1.0/rcl-rrl1+
835  (1.0-amol)*std::log(rcl*rrl1))/EE;
836  }
837 
838  //Total cross section per molecule for the active shell, in cm2
839  G4double XHTOT = XHC + XHDL + XHDT;
840 
841  //very small cross section, do nothing
842  if (XHTOT < 1.e-14*barn)
843  {
844  kineticEnergy1=kineticEnergy;
845  cosThetaPrimary=1.0;
846  energySecondary=0.0;
847  cosThetaSecondary=1.0;
848  targetOscillator = numberOfOscillators-1;
849  return;
850  }
851 
852  //decide which kind of interaction we'll have
853  TST = XHTOT*G4UniformRand();
854 
855 
856  // Hard close collision
857  G4double TS1 = XHC;
858 
859  if (TST < TS1)
860  {
861  G4double A = 5.0*amol;
862  G4double ARCL = A*0.5*rcl;
863  G4double rk=0.;
864  G4bool loopAgain = false;
865  do
866  {
867  loopAgain = false;
868  G4double fb = (1.0+ARCL)*G4UniformRand();
869  if (fb < 1)
870  rk = rcl/(1.0-fb*(1.0-(rcl+rcl)));
871  else
872  rk = rcl + (fb-1.0)*(0.5-rcl)/ARCL;
873  G4double rk2 = rk*rk;
874  G4double rkf = rk/(1.0-rk);
875  G4double phi = 1.0+rkf*rkf-rkf+amol*(rk2+rkf);
876  if (G4UniformRand()*(1.0+A*rk2) > phi)
877  loopAgain = true;
878  }while(loopAgain);
879  //energy and scattering angle (primary electron)
880  G4double deltaE = rk*EE;
881 
882  kineticEnergy1 = kineticEnergy - deltaE;
883  cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
884  //energy and scattering angle of the delta ray
885  energySecondary = deltaE - ionEne; //subtract ionisation energy
886  cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
887  if (verboseLevel > 3)
888  G4cout << "SampleFinalStateElectron: sampled close collision " << G4endl;
889  return;
890  }
891 
892  //Hard distant longitudinal collisions
893  TS1 += XHDL;
894  G4double deltaE = resEne;
895  kineticEnergy1 = kineticEnergy - deltaE;
896 
897  if (TST < TS1)
898  {
899  G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
900  G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
901  - (QS*0.5/electron_mass_c2));
902  G4double QTREV = Q*(Q+2.0*electron_mass_c2);
904  cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
905  if (cosThetaPrimary > 1.)
906  cosThetaPrimary = 1.0;
907  //energy and emission angle of the delta ray
908  energySecondary = deltaE - ionEne;
909  cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
910  if (cosThetaSecondary > 1.0)
911  cosThetaSecondary = 1.0;
912  if (verboseLevel > 3)
913  G4cout << "SampleFinalStateElectron: sampled distant longitudinal collision " << G4endl;
914  return;
915  }
916 
917  //Hard distant transverse collisions
918  cosThetaPrimary = 1.0;
919  //energy and emission angle of the delta ray
920  energySecondary = deltaE - ionEne;
921  cosThetaSecondary = 0.5;
922  if (verboseLevel > 3)
923  G4cout << "SampleFinalStateElectron: sampled distant transverse collision " << G4endl;
924 
925  return;
926 }
927 
928 
929 
930 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
932  G4double cutEnergy,
933  G4double kineticEnergy)
934 {
935  // This method sets the final ionisation parameters
936  // kineticEnergy1, cosThetaPrimary (= updates of the primary e-)
937  // energySecondary, cosThetaSecondary (= info of delta-ray)
938  // targetOscillator (= ionised oscillator)
939  //
940  // The method implements SUBROUTINE PINa of Penelope
941  //
942 
944  size_t numberOfOscillators = theTable->size();
945  const G4PenelopeCrossSection* theXS =
947  cutEnergy);
949 
950  // Selection of the active oscillator
951  G4double TST = G4UniformRand();
952  targetOscillator = numberOfOscillators-1; //initialization, last oscillator
953  G4double XSsum = 0.;
954  for (size_t i=0;i<numberOfOscillators-1;i++)
955  {
956  XSsum += theXS->GetNormalizedShellCrossSection(i,kineticEnergy);
957  if (XSsum > TST)
958  {
959  targetOscillator = (G4int) i;
960  break;
961  }
962  }
963 
964  if (verboseLevel > 3)
965  {
966  G4cout << "SampleFinalStatePositron: sampled oscillator #" << targetOscillator << "." << G4endl;
967  G4cout << "Ionisation energy: " << (*theTable)[targetOscillator]->GetIonisationEnergy()/eV
968  << " eV " << G4endl;
969  G4cout << "Resonance energy: : " << (*theTable)[targetOscillator]->GetResonanceEnergy()/eV
970  << " eV " << G4endl;
971  }
972 
973  //Constants
974  G4double rb = kineticEnergy + 2.0*electron_mass_c2;
975  G4double gam = 1.0+kineticEnergy/electron_mass_c2;
976  G4double gam2 = gam*gam;
977  G4double beta2 = (gam2-1.0)/gam2;
978  G4double g12 = (gam+1.0)*(gam+1.0);
979  G4double amol = ((gam-1.0)/gam)*((gam-1.0)/gam);
980  //Bhabha coefficients
981  G4double bha1 = amol*(2.0*g12-1.0)/(gam2-1.0);
982  G4double bha2 = amol*(3.0+1.0/g12);
983  G4double bha3 = amol*2.0*gam*(gam-1.0)/g12;
984  G4double bha4 = amol*(gam-1.0)*(gam-1.0)/g12;
985 
986  //
987  //Partial cross section of the active oscillator
988  //
989  G4double resEne = (*theTable)[targetOscillator]->GetResonanceEnergy();
990  G4double invResEne = 1.0/resEne;
991  G4double ionEne = (*theTable)[targetOscillator]->GetIonisationEnergy();
992  G4double cutoffEne = (*theTable)[targetOscillator]->GetCutoffRecoilResonantEnergy();
993 
994  G4double XHDL = 0.;
995  G4double XHDT = 0.;
996  G4double QM = 0.;
997  G4double cps = 0.;
998  G4double cp = 0.;
999 
1000  //Distant excitations XS (same as for electrons)
1001  if (resEne > cutEnergy && resEne < kineticEnergy)
1002  {
1003  cps = kineticEnergy*rb;
1004  cp = std::sqrt(cps);
1005  G4double XHDT0 = std::max(std::log(gam2)-beta2-delta,0.);
1006  if (resEne > 1.0e-6*kineticEnergy)
1007  {
1008  G4double cpp = std::sqrt((kineticEnergy-resEne)*(kineticEnergy-resEne+2.0*electron_mass_c2));
1009  QM = std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2;
1010  }
1011  else
1012  {
1013  QM = resEne*resEne/(beta2*2.0*electron_mass_c2);
1014  QM *= (1.0-QM*0.5/electron_mass_c2);
1015  }
1016  if (QM < cutoffEne)
1017  {
1018  XHDL = std::log(cutoffEne*(QM+2.0*electron_mass_c2)/(QM*(cutoffEne+2.0*electron_mass_c2)))
1019  *invResEne;
1020  XHDT = XHDT0*invResEne;
1021  }
1022  else
1023  {
1024  QM = cutoffEne;
1025  XHDL = 0.;
1026  XHDT = 0.;
1027  }
1028  }
1029  else
1030  {
1031  QM = cutoffEne;
1032  cps = 0.;
1033  cp = 0.;
1034  XHDL = 0.;
1035  XHDT = 0.;
1036  }
1037  //Close collisions (Bhabha)
1038  G4double wmaxc = kineticEnergy;
1039  G4double wcl = std::max(cutEnergy,cutoffEne);
1040  G4double rcl = wcl/kineticEnergy;
1041  G4double XHC = 0.;
1042  if (wcl < wmaxc)
1043  {
1044  G4double rl1 = 1.0-rcl;
1045  XHC = ((1.0/rcl-1.0)+bha1*std::log(rcl)+bha2*rl1
1046  + (bha3/2.0)*(rcl*rcl-1.0)
1047  + (bha4/3.0)*(1.0-rcl*rcl*rcl))/kineticEnergy;
1048  }
1049 
1050  //Total cross section per molecule for the active shell, in cm2
1051  G4double XHTOT = XHC + XHDL + XHDT;
1052 
1053  //very small cross section, do nothing
1054  if (XHTOT < 1.e-14*barn)
1055  {
1056  kineticEnergy1=kineticEnergy;
1057  cosThetaPrimary=1.0;
1058  energySecondary=0.0;
1059  cosThetaSecondary=1.0;
1060  targetOscillator = numberOfOscillators-1;
1061  return;
1062  }
1063 
1064  //decide which kind of interaction we'll have
1065  TST = XHTOT*G4UniformRand();
1066 
1067  // Hard close collision
1068  G4double TS1 = XHC;
1069  if (TST < TS1)
1070  {
1071  G4double rl1 = 1.0-rcl;
1072  G4double rk=0.;
1073  G4bool loopAgain = false;
1074  do
1075  {
1076  loopAgain = false;
1077  rk = rcl/(1.0-G4UniformRand()*rl1);
1078  G4double phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1079  if (G4UniformRand() > phi)
1080  loopAgain = true;
1081  }while(loopAgain);
1082  //energy and scattering angle (primary electron)
1083  G4double deltaE = rk*kineticEnergy;
1084  kineticEnergy1 = kineticEnergy - deltaE;
1085  cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kineticEnergy*(rb-deltaE)));
1086  //energy and scattering angle of the delta ray
1087  energySecondary = deltaE - ionEne; //subtract ionisation energy
1088  cosThetaSecondary= std::sqrt(deltaE*rb/(kineticEnergy*(deltaE+2.0*electron_mass_c2)));
1089  if (verboseLevel > 3)
1090  G4cout << "SampleFinalStatePositron: sampled close collision " << G4endl;
1091  return;
1092  }
1093 
1094  //Hard distant longitudinal collisions
1095  TS1 += XHDL;
1096  G4double deltaE = resEne;
1097  kineticEnergy1 = kineticEnergy - deltaE;
1098  if (TST < TS1)
1099  {
1100  G4double QS = QM/(1.0+QM*0.5/electron_mass_c2);
1101  G4double Q = QS/(std::pow((QS/cutoffEne)*(1.0+cutoffEne*0.5/electron_mass_c2),G4UniformRand())
1102  - (QS*0.5/electron_mass_c2));
1103  G4double QTREV = Q*(Q+2.0*electron_mass_c2);
1105  cosThetaPrimary = (cpps+cps-QTREV)/(2.0*cp*std::sqrt(cpps));
1106  if (cosThetaPrimary > 1.)
1107  cosThetaPrimary = 1.0;
1108  //energy and emission angle of the delta ray
1109  energySecondary = deltaE - ionEne;
1110  cosThetaSecondary = 0.5*(deltaE*(kineticEnergy+rb-deltaE)+QTREV)/std::sqrt(cps*QTREV);
1111  if (cosThetaSecondary > 1.0)
1112  cosThetaSecondary = 1.0;
1113  if (verboseLevel > 3)
1114  G4cout << "SampleFinalStatePositron: sampled distant longitudinal collision " << G4endl;
1115  return;
1116  }
1117 
1118  //Hard distant transverse collisions
1119  cosThetaPrimary = 1.0;
1120  //energy and emission angle of the delta ray
1121  energySecondary = deltaE - ionEne;
1122  cosThetaSecondary = 0.5;
1123 
1124  if (verboseLevel > 3)
1125  G4cout << "SampleFinalStatePositron: sampled distant transverse collision " << G4endl;
1126 
1127  return;
1128 }
1129 
1130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1131 
1133 {
1134  if(!fParticle) {
1135  fParticle = p;
1136  }
1137 }