Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups 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$
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 
54 #include "G4LivermoreIonisationModel.hh"
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"
62 #include "G4ParticleChangeForLoss.hh"
63 #include "G4Electron.hh"
64 #include "G4CrossSectionHandler.hh"
65 #include "G4VEMDataSet.hh"
66 #include "G4eIonisationCrossSectionHandler.hh"
67 #include "G4eIonisationSpectrum.hh"
68 #include "G4VEnergySpectrum.hh"
69 #include "G4SemiLogInterpolation.hh"
70 #include "G4AtomicTransitionManager.hh"
71 #include "G4AtomicShell.hh"
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
75 
76 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
77  const G4String& nam) :
78  G4VEmModel(nam), fParticleChange(0),
79  isInitialised(false),
80  crossSectionHandler(0), energySpectrum(0)
81 {
82  fIntrinsicLowEnergyLimit = 10.0*eV;
83  fIntrinsicHighEnergyLimit = 100.0*GeV;
84 
85  verboseLevel = 0;
86 
87  transitionManager = G4AtomicTransitionManager::Instance();
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
92 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
93 {
94  delete energySpectrum;
95  delete crossSectionHandler;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
99 
100 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
101  const G4DataVector& cuts)
102 {
103  //Check that the Livermore Ionisation is NOT attached to e+
104  if (particle != G4Electron::Electron())
105  {
106  G4Exception("G4LivermoreIonisationModel::Initialise",
107  "em0002",FatalException,"Livermore Ionisation Model is applicable only to electrons");
108  }
109 
110  //Read energy spectrum
111  if (energySpectrum)
112  {
113  delete energySpectrum;
114  energySpectrum = 0;
115  }
116  energySpectrum = new G4eIonisationSpectrum();
117  if (verboseLevel > 3)
118  G4cout << "G4VEnergySpectrum is initialized" << G4endl;
119 
120  //Initialize cross section handler
121  if (crossSectionHandler)
122  {
123  delete crossSectionHandler;
124  crossSectionHandler = 0;
125  }
126 
127  const size_t nbins = 20;
128  G4double emin = LowEnergyLimit();
129  G4double emax = HighEnergyLimit();
130  G4int ndec = G4int(std::log10(emax/emin) + 0.5);
131  if(ndec <= 0) { ndec = 1; }
132 
133  G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
134  crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
135  emin,emax,nbins*ndec);
136  crossSectionHandler->Clear();
137  crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
138  //This is used to retrieve cross section values later on
139  G4VEMDataSet* emdata =
140  crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
141  //The method BuildMeanFreePathForMaterials() is required here only to force
142  //the building of an internal table: the output pointer can be deleted
143  delete emdata;
144 
145  if (verboseLevel > 0)
146  {
147  G4cout << "Livermore Ionisation model is initialized " << G4endl
148  << "Energy range: "
149  << LowEnergyLimit() / keV << " keV - "
150  << HighEnergyLimit() / GeV << " GeV"
151  << G4endl;
152  }
153 
154  if (verboseLevel > 3)
155  {
156  G4cout << "Cross section data: " << G4endl;
157  crossSectionHandler->PrintData();
158  G4cout << "Parameters: " << G4endl;
159  energySpectrum->PrintData();
160  }
161 
162  if(isInitialised) { return; }
163  fParticleChange = GetParticleChangeForLoss();
164  isInitialised = true;
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168 
169 G4double
170 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
171  G4double energy,
172  G4double Z, G4double,
173  G4double cutEnergy,
174  G4double)
175 {
176  G4int iZ = (G4int) Z;
177  if (!crossSectionHandler)
178  {
179  G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
180  "em1007",FatalException,"The cross section handler is not correctly initialized");
181  return 0;
182  }
183 
184  //The cut is already included in the crossSectionHandler
185  G4double cs =
186  crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
187 
188  if (verboseLevel > 1)
189  {
190  G4cout << "G4LivermoreIonisationModel " << G4endl;
191  G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
192  energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
193  }
194  return cs;
195 }
196 
197 
198 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
199 
200 G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
201  const G4ParticleDefinition* ,
202  G4double kineticEnergy,
203  G4double cutEnergy)
204 {
205  G4double sPower = 0.0;
206 
207  const G4ElementVector* theElementVector = material->GetElementVector();
208  size_t NumberOfElements = material->GetNumberOfElements() ;
209  const G4double* theAtomicNumDensityVector =
210  material->GetAtomicNumDensityVector();
211 
212  // loop for elements in the material
213  for (size_t iel=0; iel<NumberOfElements; iel++ )
214  {
215  G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
216  G4int nShells = transitionManager->NumberOfShells(iZ);
217  for (G4int n=0; n<nShells; n++)
218  {
219  G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
220  kineticEnergy, n);
221  G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
222  sPower += e * cs * theAtomicNumDensityVector[iel];
223  }
224  G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
225  sPower += esp * theAtomicNumDensityVector[iel];
226  }
227 
228  if (verboseLevel > 2)
229  {
230  G4cout << "G4LivermoreIonisationModel " << G4endl;
231  G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
232  kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
233  }
234 
235  return sPower;
236 }
237 
238 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
239 
240 void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
241  const G4MaterialCutsCouple* couple,
242  const G4DynamicParticle* aDynamicParticle,
243  G4double cutE,
244  G4double maxE)
245 {
246 
247  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
248 
249  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
250  {
251  fParticleChange->SetProposedKineticEnergy(0.);
252  fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
253  return;
254  }
255 
256  // Select atom and shell
257  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
258  G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
259  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
260  G4double bindingEnergy = shell->BindingEnergy();
261 
262  // Sample delta energy using energy interval for delta-electrons
263  G4double energyMax =
264  std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
265  G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
266  kineticEnergy, shellIndex);
267 
268  if (energyDelta == 0.) //nothing happens
269  { return; }
270 
271  // Transform to shell potential
272  G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
273  G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
274 
275  // sampling of scattering angle neglecting atomic motion
276  G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
277  G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
278 
279  G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
280  / (deltaMom * primaryMom);
281  if (cost > 1.) { cost = 1.; }
282  G4double sint = std::sqrt((1. - cost)*(1. + cost));
283  G4double phi = twopi * G4UniformRand();
284  G4double dirx = sint * std::cos(phi);
285  G4double diry = sint * std::sin(phi);
286  G4double dirz = cost;
287 
288  // Rotate to incident electron direction
289  G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
290  G4ThreeVector deltaDir(dirx,diry,dirz);
291  deltaDir.rotateUz(primaryDirection);
292  //Updated components
293  dirx = deltaDir.x();
294  diry = deltaDir.y();
295  dirz = deltaDir.z();
296 
297  // Take into account atomic motion del is relative momentum of the motion
298  // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
299  cost = 2.0*G4UniformRand() - 1.0;
300  sint = std::sqrt(1. - cost*cost);
301  phi = twopi * G4UniformRand();
302  G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
303  / deltaMom;
304  dirx += del* sint * std::cos(phi);
305  diry += del* sint * std::sin(phi);
306  dirz += del* cost;
307 
308  // Find out new primary electron direction
309  G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
310  G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
311  G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
312 
313  //Ok, ready to create the delta ray
314  G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
315  theDeltaRay->SetKineticEnergy(energyDelta);
316  G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
317  dirx *= norm;
318  diry *= norm;
319  dirz *= norm;
320  theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
321  theDeltaRay->SetDefinition(G4Electron::Electron());
322  fvect->push_back(theDeltaRay);
323 
324  //This is the amount of energy available for fluorescence
325  G4double theEnergyDeposit = bindingEnergy;
326 
327  // fill ParticleChange
328  // changed energy and momentum of the actual particle
329  G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
330  if(finalKinEnergy < 0.0)
331  {
332  theEnergyDeposit += finalKinEnergy;
333  finalKinEnergy = 0.0;
334  }
335  else
336  {
337  G4double normLocal = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
338  finalPx *= normLocal;
339  finalPy *= normLocal;
340  finalPz *= normLocal;
341  fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
342  }
343  fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
344 
345  if (theEnergyDeposit < 0)
346  {
347  G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
348  << theEnergyDeposit/eV << " eV" << G4endl;
349  theEnergyDeposit = 0.0;
350  }
351 
352  //Assign local energy deposit
353  fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
354 
355  if (verboseLevel > 1)
356  {
357  G4cout << "-----------------------------------------------------------" << G4endl;
358  G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
359  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
360  G4cout << "-----------------------------------------------------------" << G4endl;
361  G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
362  G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
363  G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
364  G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
365  G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
366  << " keV" << G4endl;
367  G4cout << "-----------------------------------------------------------" << G4endl;
368  }
369  return;
370 }
371 
372 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373 
Generated on Sat May 25 2013 14:33:34 for Geant4 by   doxygen 1.8.4