Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4LivermoreIonisationModel.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: G4LivermoreIonisationModel.cc 93810 2015-11-02 11:27:56Z gcosmo $
27 //
28 // Author: Luciano Pandola
29 // on base of G4LowEnergyIonisation developed by A.Forti and V.Ivanchenko
30 //
31 // History:
32 // --------
33 // 12 Jan 2009 L Pandola Migration from process to model
34 // 03 Mar 2009 L Pandola Bug fix (release memory in the destructor)
35 // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
36 // - apply internal high-energy limit only in constructor
37 // - do not apply low-energy limit (default is 0)
38 // - simplify sampling of deexcitation by using cut in energy
39 // - set activation of Auger "false"
40 // - remove initialisation of element selectors
41 // 19 May 2009 L Pandola Explicitely set to zero pointers deleted in
42 // Initialise(), since they might be checked later on
43 // 23 Oct 2009 L Pandola
44 // - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
45 // set as "true" (default would be false)
46 // 12 Oct 2010 L Pandola
47 // - add debugging information about energy in
48 // SampleDeexcitationAlongStep()
49 // - generate fluorescence SampleDeexcitationAlongStep() only above
50 // the cuts.
51 // 01 Jun 2011 V Ivanchenko general cleanup - all old deexcitation code removed
52 //
53 
55 #include "G4PhysicalConstants.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4ParticleDefinition.hh"
58 #include "G4MaterialCutsCouple.hh"
59 #include "G4ProductionCutsTable.hh"
60 #include "G4DynamicParticle.hh"
61 #include "G4Element.hh"
63 #include "G4Electron.hh"
64 #include "G4CrossSectionHandler.hh"
65 #include "G4VEMDataSet.hh"
67 #include "G4eIonisationSpectrum.hh"
68 #include "G4VEnergySpectrum.hh"
71 #include "G4AtomicShell.hh"
72 #include "G4DeltaAngle.hh"
73 
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 
76 
78  const G4String& nam) :
79  G4VEmModel(nam), fParticleChange(0),
80  isInitialised(false),
81  crossSectionHandler(0), energySpectrum(0)
82 {
83  fIntrinsicLowEnergyLimit = 12.*eV;
84  fIntrinsicHighEnergyLimit = 100.0*GeV;
85 
86  verboseLevel = 0;
88 
89  transitionManager = G4AtomicTransitionManager::Instance();
90 }
91 
92 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
93 
95 {
96  delete energySpectrum;
97  delete crossSectionHandler;
98 }
99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
101 
103  const G4DataVector& cuts)
104 {
105  //Check that the Livermore Ionisation is NOT attached to e+
106  if (particle != G4Electron::Electron())
107  {
108  G4Exception("G4LivermoreIonisationModel::Initialise",
109  "em0002",FatalException,
110  "Livermore Ionisation Model is applicable only to electrons");
111  }
112 
113  transitionManager->Initialise();
114 
115  //Read energy spectrum
116  if (energySpectrum)
117  {
118  delete energySpectrum;
119  energySpectrum = 0;
120  }
121  energySpectrum = new G4eIonisationSpectrum();
122  if (verboseLevel > 3)
123  G4cout << "G4VEnergySpectrum is initialized" << G4endl;
124 
125  //Initialize cross section handler
126  if (crossSectionHandler)
127  {
128  delete crossSectionHandler;
129  crossSectionHandler = 0;
130  }
131 
132  const size_t nbins = 20;
133  G4double emin = LowEnergyLimit();
135  G4int ndec = G4int(std::log10(emax/emin) + 0.5);
136  if(ndec <= 0) { ndec = 1; }
137 
138  G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
139  crossSectionHandler =
140  new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
141  emin,emax,nbins*ndec);
142  crossSectionHandler->Clear();
143  crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
144  //This is used to retrieve cross section values later on
145  G4VEMDataSet* emdata =
146  crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
147  //The method BuildMeanFreePathForMaterials() is required here only to force
148  //the building of an internal table: the output pointer can be deleted
149  delete emdata;
150 
151  if (verboseLevel > 0)
152  {
153  G4cout << "Livermore Ionisation model is initialized " << G4endl
154  << "Energy range: "
155  << LowEnergyLimit() / keV << " keV - "
156  << HighEnergyLimit() / GeV << " GeV"
157  << G4endl;
158  }
159 
160  if (verboseLevel > 3)
161  {
162  G4cout << "Cross section data: " << G4endl;
163  crossSectionHandler->PrintData();
164  G4cout << "Parameters: " << G4endl;
165  energySpectrum->PrintData();
166  }
167 
168  if(isInitialised) { return; }
170  isInitialised = true;
171 }
172 
173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174 
175 G4double
177  const G4ParticleDefinition*,
180  G4double cutEnergy,
181  G4double)
182 {
183  G4int iZ = (G4int) Z;
184  if (!crossSectionHandler)
185  {
186  G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
187  "em1007",FatalException,
188  "The cross section handler is not correctly initialized");
189  return 0;
190  }
191 
192  //The cut is already included in the crossSectionHandler
193  G4double cs =
194  crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,
195  cutEnergy,
196  iZ);
197 
198  if (verboseLevel > 1)
199  {
200  G4cout << "G4LivermoreIonisationModel " << G4endl;
201  G4cout << "Cross section for delta emission > "
202  << cutEnergy/keV << " keV at "
203  << energy/keV << " keV and Z = " << iZ << " --> "
204  << cs/barn << " barn" << G4endl;
205  }
206  return cs;
207 }
208 
209 
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211 
212 G4double
214  const G4ParticleDefinition*,
215  G4double kineticEnergy,
216  G4double cutEnergy)
217 {
218  G4double sPower = 0.0;
219 
220  const G4ElementVector* theElementVector = material->GetElementVector();
221  size_t NumberOfElements = material->GetNumberOfElements() ;
222  const G4double* theAtomicNumDensityVector =
223  material->GetAtomicNumDensityVector();
224 
225  // loop for elements in the material
226  for (size_t iel=0; iel<NumberOfElements; iel++ )
227  {
228  G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
229  G4int nShells = transitionManager->NumberOfShells(iZ);
230  for (G4int n=0; n<nShells; n++)
231  {
232  G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
233  kineticEnergy, n);
234  G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
235  sPower += e * cs * theAtomicNumDensityVector[iel];
236  }
237  G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
238  sPower += esp * theAtomicNumDensityVector[iel];
239  }
240 
241  if (verboseLevel > 2)
242  {
243  G4cout << "G4LivermoreIonisationModel " << G4endl;
244  G4cout << "Stopping power < " << cutEnergy/keV
245  << " keV at " << kineticEnergy/keV << " keV = "
246  << sPower/(keV/mm) << " keV/mm" << G4endl;
247  }
248 
249  return sPower;
250 }
251 
252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
253 
255  std::vector<G4DynamicParticle*>* fvect,
256  const G4MaterialCutsCouple* couple,
257  const G4DynamicParticle* aDynamicParticle,
258  G4double cutE,
259  G4double maxE)
260 {
261 
262  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
263 
264  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
265  {
268  return;
269  }
270 
271  // Select atom and shell
272  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
273  G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
274  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
276 
277  // Sample delta energy using energy interval for delta-electrons
278  G4double energyMax =
279  std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
280  G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
281  kineticEnergy, shellIndex);
282 
283  if (energyDelta == 0.) //nothing happens
284  { return; }
285 
287  G4DynamicParticle* delta = new G4DynamicParticle(electron,
288  GetAngularDistribution()->SampleDirectionForShell(aDynamicParticle, energyDelta,
289  Z, shellIndex,
290  couple->GetMaterial()),
291  energyDelta);
292 
293  fvect->push_back(delta);
294 
295  // Change kinematics of primary particle
296  G4ThreeVector direction = aDynamicParticle->GetMomentumDirection();
297  G4double totalMomentum = std::sqrt(kineticEnergy*(kineticEnergy + 2*electron_mass_c2));
298 
299  G4ThreeVector finalP = totalMomentum*direction - delta->GetMomentum();
300  finalP = finalP.unit();
301 
302  //This is the amount of energy available for fluorescence
303  G4double theEnergyDeposit = bindingEnergy;
304 
305  // fill ParticleChange
306  // changed energy and momentum of the actual particle
307  G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
308  if(finalKinEnergy < 0.0)
309  {
310  theEnergyDeposit += finalKinEnergy;
311  finalKinEnergy = 0.0;
312  }
313  else
314  {
316  }
317  fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
318 
319  if (theEnergyDeposit < 0)
320  {
321  G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
322  << theEnergyDeposit/eV << " eV" << G4endl;
323  theEnergyDeposit = 0.0;
324  }
325 
326  //Assign local energy deposit
327  fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
328 
329  if (verboseLevel > 1)
330  {
331  G4cout << "-----------------------------------------------------------" << G4endl;
332  G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
333  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
334  G4cout << "-----------------------------------------------------------" << G4endl;
335  G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
336  G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
337  G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
338  G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
339  G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
340  << " keV" << G4endl;
341  G4cout << "-----------------------------------------------------------" << G4endl;
342  }
343  return;
344 }
345 
346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
347 
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:640
virtual void PrintData() const =0
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:118
static constexpr double mm
Definition: G4SIunits.hh:115
std::vector< G4Element * > G4ElementVector
G4double GetCrossSectionAboveThresholdForElement(G4double energy, G4double cutEnergy, G4int Z)
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:633
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4int SelectRandomShell(G4int Z, G4double e) const
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:616
G4double BindingEnergy() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
static constexpr double electron_mass_c2
G4ParticleChangeForLoss * fParticleChange
virtual G4double SampleEnergy(G4int Z, G4double minKineticEnergy, G4double maxKineticEnergy, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const =0
G4double FindValue(G4int Z, G4double e) const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
virtual G4double AverageEnergy(G4int Z, G4double minKineticEnergy, G4double maxKineticEnergy, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const =0
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static constexpr double eV
Definition: G4SIunits.hh:215
G4VEMDataSet * BuildMeanFreePathForMaterials(const G4DataVector *energyCuts=0)
virtual G4double MaxEnergyOfSecondaries(G4double kineticEnergy, G4int Z=0, const G4ParticleDefinition *pd=0) const =0
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:216
virtual G4double ComputeDEDXPerVolume(const G4Material *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy)
void LoadShellData(const G4String &dataFile)
G4LivermoreIonisationModel(const G4ParticleDefinition *p=0, const G4String &processName="LowEnergyIoni")
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static const G4double emax
G4double energy(const ThreeVector &p, const G4double m)
Hep3Vector unit() const
virtual G4double Excitation(G4int Z, G4double kineticEnergy) const =0
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:623
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double GeV
Definition: G4SIunits.hh:217
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
double G4double
Definition: G4Types.hh:76
static constexpr double barn
Definition: G4SIunits.hh:105
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
G4int SelectRandomAtom(const G4MaterialCutsCouple *couple, G4double e) const
G4ThreeVector GetMomentum() const
const G4Material * GetMaterial() const