Geant4_10
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 74822 2013-10-22 14:42:13Z 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 = 10.0*eV;
84  fIntrinsicHighEnergyLimit = 100.0*GeV;
85 
86  verboseLevel = 0;
87  //SetAngularDistribution(new G4DeltaAngle());
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  //Read energy spectrum
114  if (energySpectrum)
115  {
116  delete energySpectrum;
117  energySpectrum = 0;
118  }
119  energySpectrum = new G4eIonisationSpectrum();
120  if (verboseLevel > 3)
121  G4cout << "G4VEnergySpectrum is initialized" << G4endl;
122 
123  //Initialize cross section handler
124  if (crossSectionHandler)
125  {
126  delete crossSectionHandler;
127  crossSectionHandler = 0;
128  }
129 
130  const size_t nbins = 20;
131  G4double emin = LowEnergyLimit();
132  G4double emax = HighEnergyLimit();
133  G4int ndec = G4int(std::log10(emax/emin) + 0.5);
134  if(ndec <= 0) { ndec = 1; }
135 
136  G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
137  crossSectionHandler =
138  new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
139  emin,emax,nbins*ndec);
140  crossSectionHandler->Clear();
141  crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
142  //This is used to retrieve cross section values later on
143  G4VEMDataSet* emdata =
144  crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
145  //The method BuildMeanFreePathForMaterials() is required here only to force
146  //the building of an internal table: the output pointer can be deleted
147  delete emdata;
148 
149  if (verboseLevel > 0)
150  {
151  G4cout << "Livermore Ionisation model is initialized " << G4endl
152  << "Energy range: "
153  << LowEnergyLimit() / keV << " keV - "
154  << HighEnergyLimit() / GeV << " GeV"
155  << G4endl;
156  }
157 
158  if (verboseLevel > 3)
159  {
160  G4cout << "Cross section data: " << G4endl;
161  crossSectionHandler->PrintData();
162  G4cout << "Parameters: " << G4endl;
163  energySpectrum->PrintData();
164  }
165 
166  if(isInitialised) { return; }
168  isInitialised = true;
169 }
170 
171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
172 
173 G4double
175  const G4ParticleDefinition*,
178  G4double cutEnergy,
179  G4double)
180 {
181  G4int iZ = (G4int) Z;
182  if (!crossSectionHandler)
183  {
184  G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
185  "em1007",FatalException,
186  "The cross section handler is not correctly initialized");
187  return 0;
188  }
189 
190  //The cut is already included in the crossSectionHandler
191  G4double cs =
192  crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,
193  cutEnergy,
194  iZ);
195 
196  if (verboseLevel > 1)
197  {
198  G4cout << "G4LivermoreIonisationModel " << G4endl;
199  G4cout << "Cross section for delta emission > "
200  << cutEnergy/keV << " keV at "
201  << energy/keV << " keV and Z = " << iZ << " --> "
202  << cs/barn << " barn" << G4endl;
203  }
204  return cs;
205 }
206 
207 
208 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
209 
210 G4double
212  const G4ParticleDefinition*,
213  G4double kineticEnergy,
214  G4double cutEnergy)
215 {
216  G4double sPower = 0.0;
217 
218  const G4ElementVector* theElementVector = material->GetElementVector();
219  size_t NumberOfElements = material->GetNumberOfElements() ;
220  const G4double* theAtomicNumDensityVector =
221  material->GetAtomicNumDensityVector();
222 
223  // loop for elements in the material
224  for (size_t iel=0; iel<NumberOfElements; iel++ )
225  {
226  G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
227  G4int nShells = transitionManager->NumberOfShells(iZ);
228  for (G4int n=0; n<nShells; n++)
229  {
230  G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
231  kineticEnergy, n);
232  G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
233  sPower += e * cs * theAtomicNumDensityVector[iel];
234  }
235  G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
236  sPower += esp * theAtomicNumDensityVector[iel];
237  }
238 
239  if (verboseLevel > 2)
240  {
241  G4cout << "G4LivermoreIonisationModel " << G4endl;
242  G4cout << "Stopping power < " << cutEnergy/keV
243  << " keV at " << kineticEnergy/keV << " keV = "
244  << sPower/(keV/mm) << " keV/mm" << G4endl;
245  }
246 
247  return sPower;
248 }
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251 
253  std::vector<G4DynamicParticle*>* fvect,
254  const G4MaterialCutsCouple* couple,
255  const G4DynamicParticle* aDynamicParticle,
256  G4double cutE,
257  G4double maxE)
258 {
259 
260  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
261 
262  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
263  {
266  return;
267  }
268 
269  // Select atom and shell
270  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
271  G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
272  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
274 
275  // Sample delta energy using energy interval for delta-electrons
276  G4double energyMax =
277  std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
278  G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
279  kineticEnergy, shellIndex);
280 
281  if (energyDelta == 0.) //nothing happens
282  { return; }
283 
284  // Transform to shell potential
285  G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
286  G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
287 
288  // sampling of scattering angle neglecting atomic motion
289  G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
290  G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
291 
292  G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
293  / (deltaMom * primaryMom);
294  if (cost > 1.) { cost = 1.; }
295  G4double sint = std::sqrt((1. - cost)*(1. + cost));
296  G4double phi = twopi * G4UniformRand();
297  G4double dirx = sint * std::cos(phi);
298  G4double diry = sint * std::sin(phi);
299  G4double dirz = cost;
300 
301  // Rotate to incident electron direction
302  G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
303  G4ThreeVector deltaDir(dirx,diry,dirz);
304  deltaDir.rotateUz(primaryDirection);
305  //Updated components
306  dirx = deltaDir.x();
307  diry = deltaDir.y();
308  dirz = deltaDir.z();
309 
310  // Take into account atomic motion del is relative momentum of the motion
311  // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
312  cost = 2.0*G4UniformRand() - 1.0;
313  sint = std::sqrt(1. - cost*cost);
314  phi = twopi * G4UniformRand();
315  G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
316  / deltaMom;
317  dirx += del* sint * std::cos(phi);
318  diry += del* sint * std::sin(phi);
319  dirz += del* cost;
320 
321  // Find out new primary electron direction
322  G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
323  G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
324  G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
325 
326  //Ok, ready to create the delta ray
327  G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
328  theDeltaRay->SetKineticEnergy(energyDelta);
329  G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
330  dirx *= norm;
331  diry *= norm;
332  dirz *= norm;
333  theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
334  theDeltaRay->SetDefinition(G4Electron::Electron());
335  fvect->push_back(theDeltaRay);
336 
337  //This is the amount of energy available for fluorescence
338  G4double theEnergyDeposit = bindingEnergy;
339 
340  // fill ParticleChange
341  // changed energy and momentum of the actual particle
342  G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
343  if(finalKinEnergy < 0.0)
344  {
345  theEnergyDeposit += finalKinEnergy;
346  finalKinEnergy = 0.0;
347  }
348  else
349  {
350  G4double normLocal = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
351  finalPx *= normLocal;
352  finalPy *= normLocal;
353  finalPz *= normLocal;
354  fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
355  }
356  fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
357 
358  if (theEnergyDeposit < 0)
359  {
360  G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
361  << theEnergyDeposit/eV << " eV" << G4endl;
362  theEnergyDeposit = 0.0;
363  }
364 
365  //Assign local energy deposit
366  fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
367 
368  if (verboseLevel > 1)
369  {
370  G4cout << "-----------------------------------------------------------" << G4endl;
371  G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
372  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
373  G4cout << "-----------------------------------------------------------" << G4endl;
374  G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
375  G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
376  G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
377  G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
378  G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
379  << " keV" << G4endl;
380  G4cout << "-----------------------------------------------------------" << G4endl;
381  }
382  return;
383 }
384 
385 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386 
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
virtual void PrintData() const =0
G4ParticleChangeForLoss * GetParticleChangeForLoss()
Definition: G4VEmModel.cc:107
std::vector< G4Element * > G4ElementVector
G4double GetCrossSectionAboveThresholdForElement(G4double energy, G4double cutEnergy, G4int Z)
G4double GetKineticEnergy() const
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
double x() const
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
G4int SelectRandomShell(G4int Z, G4double e) const
Float_t norm
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4double BindingEnergy() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
double z() const
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4ParticleChangeForLoss * fParticleChange
virtual G4double SampleEnergy(G4int Z, G4double minKineticEnergy, G4double maxKineticEnergy, G4double kineticEnergy, G4int shell=0, const G4ParticleDefinition *pd=0) const =0
string material
Definition: eplot.py:19
G4double FindValue(G4int Z, G4double e) const
Char_t n[5]
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
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
Float_t Z
Definition: plot.C:39
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
void SetProposedKineticEnergy(G4double proposedKinEnergy)
float electron_mass_c2
Definition: hepunit.py:274
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:214
void SetKineticEnergy(G4double aEnergy)
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
virtual G4double Excitation(G4int Z, G4double kineticEnergy) const =0
double y() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
double G4double
Definition: G4Types.hh:76
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double bindingEnergy(G4int A, G4int Z)
static G4AtomicTransitionManager * Instance()
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
G4int SelectRandomAtom(const G4MaterialCutsCouple *couple, G4double e) const