source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeIonisationModel.cc@ 1014

Last change on this file since 1014 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 60.6 KB
RevLine 
[968]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//
[1007]26// $Id: G4PenelopeIonisationModel.cc,v 1.2 2008/12/05 09:15:43 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02 $
[968]28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 26 Nov 2008 L Pandola Migration from process to model
34//
35
36#include "G4PenelopeIonisationModel.hh"
37#include "G4ParticleDefinition.hh"
38#include "G4MaterialCutsCouple.hh"
39#include "G4ProductionCutsTable.hh"
40#include "G4DynamicParticle.hh"
41#include "G4Element.hh"
42#include "G4AtomicTransitionManager.hh"
43#include "G4AtomicDeexcitation.hh"
44#include "G4AtomicShell.hh"
45#include "G4Gamma.hh"
46#include "G4Electron.hh"
47#include "G4Positron.hh"
48#include "G4CrossSectionHandler.hh"
49#include "G4AtomicDeexcitation.hh"
50#include "G4ProcessManager.hh"
51#include "G4VEMDataSet.hh"
52
53//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
54
55
56G4PenelopeIonisationModel::G4PenelopeIonisationModel(const G4ParticleDefinition*,
57 const G4String& nam)
58 :G4VEmModel(nam),isInitialised(false),
59 kineticEnergy1(0.*eV),
60 cosThetaPrimary(1.0),
61 energySecondary(0.*eV),
62 cosThetaSecondary(0.0),
63 iOsc(-1),
64 crossSectionHandler(0),ionizationEnergy(0),
65 resonanceEnergy(0),occupationNumber(0),shellFlag(0),
66 theXSTable(0)
67{
68 fIntrinsicLowEnergyLimit = 100.0*eV;
69 fIntrinsicHighEnergyLimit = 100.0*GeV;
70 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
71 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
72 //
73
74 fUseAtomicDeexcitation = true;
75 verboseLevel= 0;
76
77 // Verbosity scale:
78 // 0 = nothing
79 // 1 = warning for energy non-conservation
80 // 2 = details of energy budget
81 // 3 = calculation of cross sections, file openings, sampling of atoms
82 // 4 = entering in methods
83
84 //These vectors do not change when materials or cut change.
85 //Therefore I can read it at the constructor
86 ionizationEnergy = new std::map<G4int,G4DataVector*>;
87 resonanceEnergy = new std::map<G4int,G4DataVector*>;
88 occupationNumber = new std::map<G4int,G4DataVector*>;
89 shellFlag = new std::map<G4int,G4DataVector*>;
90
91 ReadData(); //Read data from file
92
93}
94
95//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
96
97G4PenelopeIonisationModel::~G4PenelopeIonisationModel()
98{
99
100 if (crossSectionHandler) delete crossSectionHandler;
101
102 if (theXSTable)
103 {
104 for (size_t i=0; i<theXSTable->size(); i++)
105 delete (*theXSTable)[i];
106 delete theXSTable;
107 }
108
109
110 std::map <G4int,G4DataVector*>::iterator i;
111 for (i=ionizationEnergy->begin();i != ionizationEnergy->end();i++)
112 if (i->second) delete i->second;
113 for (i=resonanceEnergy->begin();i != resonanceEnergy->end();i++)
114 if (i->second) delete i->second;
115 for (i=occupationNumber->begin();i != occupationNumber->end();i++)
116 if (i->second) delete i->second;
117 for (i=shellFlag->begin();i != shellFlag->end();i++)
118 if (i->second) delete i->second;
119
120 if (ionizationEnergy)
121 delete ionizationEnergy;
122 if (resonanceEnergy)
123 delete resonanceEnergy;
124 if (occupationNumber)
125 delete occupationNumber;
126 if (shellFlag)
127 delete shellFlag;
128}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131
132void G4PenelopeIonisationModel::Initialise(const G4ParticleDefinition* particle,
133 const G4DataVector& cuts)
134{
135 if (verboseLevel > 3)
136 G4cout << "Calling G4PenelopeIonisationModel::Initialise()" << G4endl;
137
138 InitialiseElementSelectors(particle,cuts);
139
140 //Delete and re-initialize the cross section handler
141 if (crossSectionHandler)
142 {
143 crossSectionHandler->Clear();
144 delete crossSectionHandler;
145 }
146
147 if (theXSTable)
148 {
149 for (size_t i=0; i<theXSTable->size(); i++)
150 delete (*theXSTable)[i];
151 delete theXSTable;
152 }
153
154
155 if (LowEnergyLimit() < fIntrinsicLowEnergyLimit)
156 {
157 G4cout << "G4PenelopeIonisationModel: low energy limit increased from " <<
158 LowEnergyLimit()/eV << " eV to " << fIntrinsicLowEnergyLimit/eV << " eV" << G4endl;
159 SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
160 }
161 if (HighEnergyLimit() > fIntrinsicHighEnergyLimit)
162 {
163 G4cout << "G4PenelopeIonisationModel: high energy limit decreased from " <<
164 HighEnergyLimit()/GeV << " GeV to " << fIntrinsicHighEnergyLimit/GeV << " GeV" << G4endl;
165 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
166 }
167
168 crossSectionHandler = new G4CrossSectionHandler();
169 crossSectionHandler->Clear();
170 G4String crossSectionFile = "NULL";
171
172 if (particle == G4Electron::Electron())
173 crossSectionFile = "penelope/ion-cs-el-";
174 else if (particle == G4Positron::Positron())
175 crossSectionFile = "penelope/ion-cs-po-";
176 crossSectionHandler->LoadData(crossSectionFile);
177 //This is used to retrieve cross section values later on
178 crossSectionHandler->BuildMeanFreePathForMaterials();
179
180 if (verboseLevel > 2)
181 G4cout << "Loaded cross section files for PenelopeIonisationModel" << G4endl;
182
183
184 G4cout << "Penelope Ionisation model is initialized " << G4endl
185 << "Energy range: "
186 << LowEnergyLimit() / keV << " keV - "
187 << HighEnergyLimit() / GeV << " GeV"
188 << G4endl;
189
190 if(isInitialised) return;
191
192 if(pParticleChange)
193 fParticleChange = reinterpret_cast<G4ParticleChangeForLoss*>(pParticleChange);
194 else
195 fParticleChange = new G4ParticleChangeForLoss();
196 isInitialised = true;
197}
198
199//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
200
201G4double G4PenelopeIonisationModel::CrossSectionPerVolume(const G4Material* material,
202 const G4ParticleDefinition* theParticle,
203 G4double energy,
204 G4double cutEnergy,
205 G4double)
206{
207 // Penelope model to calculate the cross section for inelastic collisions above the
208 // threshold. It makes use of the Generalised Oscillator Strength (GOS) model from
209 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
210 //
211 // The total cross section (hard+soft) is read from a database file (element per
212 // element), while the ratio hard-to-total is calculated analytically by taking
213 // into account the atomic oscillators coming into the play for a given threshold.
214 // This is done by the method CalculateCrossSectionsRatio().
215 // For incident e- the maximum energy allowed for the delta-rays is energy/2.
216 // because particles are undistinghishable.
217 //
218 // The contribution is splitted in three parts: distant longitudinal collisions,
219 // distant transverse collisions and close collisions. Each term is described by
220 // its own analytical function.
221 // Fermi density correction is calculated analytically according to
222 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
223 //
224 if (verboseLevel > 3)
225 G4cout << "Calling CrossSectionPerVolume() of G4PenelopeIonisationModel" << G4endl;
226
227 SetupForMaterial(theParticle, material, energy);
228
229 if (!crossSectionHandler)
230 {
231 G4cout << "G4PenelopeIonisationModel::CrossSectionPerVolume" << G4endl;
232 G4cout << "The cross section handler is not correctly initialized" << G4endl;
233 G4Exception();
234 }
235
236 if (!theXSTable)
237 {
238 if (verboseLevel > 2)
239 {
240 G4cout << "G4PenelopeIonisationModel::CrossSectionPerVolume" << G4endl;
241 G4cout << "Going to build Cross Section table " << G4endl;
242 }
243 theXSTable = new std::vector<G4VEMDataSet*>;
244 theXSTable = BuildCrossSectionTable(theParticle);
245 }
246
247 G4double totalCross = 0.0;
248 G4double cross = 0.0;
249 const G4ElementVector* theElementVector = material->GetElementVector();
250 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
251 G4double electronVolumeDensity =
252 material->GetTotNbOfElectPerVolume(); //electron density
253
254 if (verboseLevel > 4)
255 G4cout << "Electron volume density of " << material->GetName() << ": " <<
256 electronVolumeDensity*cm3 << " electrons/cm3" << G4endl;
257
258 G4int nelm = material->GetNumberOfElements();
259 for (G4int i=0; i<nelm; i++)
260 {
261 G4int iZ = (G4int) (*theElementVector)[i]->GetZ();
262 G4double ratio = CalculateCrossSectionsRatio(energy,cutEnergy,iZ,electronVolumeDensity,
263 theParticle);
264 cross += theAtomNumDensityVector[i]*
265 crossSectionHandler->FindValue(iZ,energy)*ratio;
266 totalCross += theAtomNumDensityVector[i]*
267 crossSectionHandler->FindValue(iZ,energy);
268 }
269
270 if (verboseLevel > 2)
271 {
272 G4cout << "G4PenelopeIonisationModel " << G4endl;
273 G4cout << "Mean free path for delta emission > " << cutEnergy/keV << " keV at " <<
274 energy/keV << " keV = " << (1./cross)/mm << " mm" << G4endl;
275 G4cout << "Total free path for ionisation (no threshold) at " <<
276 energy/keV << " keV = " << (1./totalCross)/mm << " mm" << G4endl;
277 }
278 return cross;
279}
280
281//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282
283G4double G4PenelopeIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
284 const G4ParticleDefinition* theParticle,
285 G4double kineticEnergy,
286 G4double cutEnergy)
287{
288 // Penelope model to calculate the stopping power for soft inelastic collisions
289 // below the threshold. It makes use of the Generalised Oscillator Strength (GOS)
290 // model from
291 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
292 //
293 // The stopping power is calculated analytically using the dsigma/dW cross
294 // section from the GOS models, which includes separate contributions from
295 // distant longitudinal collisions, distant transverse collisions and
296 // close collisions. Only the atomic oscillators that come in the play
297 // (according to the threshold) are considered for the calculation.
298 // Differential cross sections have a different form for e+ and e-.
299 //
300 // Fermi density correction is calculated analytically according to
301 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
302
303 if (verboseLevel > 3)
304 G4cout << "Calling ComputeDEDX() of G4PenelopeIonisationModel" << G4endl;
305
306 G4double sPower = 0.0;
307 const G4ElementVector* theElementVector = material->GetElementVector();
308 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
309 G4double electronVolumeDensity =
310 material->GetTotNbOfElectPerVolume(); //electron density
311
312 //Loop on the elements in the material
313 G4int nelm = material->GetNumberOfElements();
314 for (G4int i=0; i<nelm; i++)
315 {
316 G4int iZ = (G4int) (*theElementVector)[i]->GetZ();
317
318 //Calculate stopping power contribution from each element
319 //Constants
320 G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
321 G4double gamma2 = gamma*gamma;
322 G4double beta2 = (gamma2-1.0)/gamma2;
323 G4double constant = pi*classic_electr_radius*classic_electr_radius
324 *2.0*electron_mass_c2/beta2;
325
326 G4double delta = CalculateDeltaFermi(kineticEnergy,iZ,electronVolumeDensity);
327 G4int nbOsc = (G4int) resonanceEnergy->find(iZ)->second->size();
328 G4double stoppingPowerForElement = 0.0;
329 //Loop on oscillators of element Z
330 for (G4int iosc=0;iosc<nbOsc;iosc++)
331 {
332 G4double S1 = 0.0;
333 G4double resEnergy = (*(resonanceEnergy->find(iZ)->second))[iosc];
334 if (theParticle == G4Electron::Electron())
335 S1 = ComputeStoppingPowerForElectrons(kineticEnergy,cutEnergy,delta,resEnergy);
336 else if (theParticle == G4Positron::Positron())
337 S1 = ComputeStoppingPowerForPositrons(kineticEnergy,cutEnergy,delta,resEnergy);
338 G4double occupNb = (*(occupationNumber->find(iZ)->second))[iosc];
339 stoppingPowerForElement += occupNb*constant*S1;
340 }
341
342 sPower += stoppingPowerForElement*theAtomNumDensityVector[i];
343 }
344 if (verboseLevel > 2)
345 {
346 G4cout << "G4PenelopeIonisationModel " << G4endl;
347 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
348 kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
349 }
350 return sPower;
351}
352
353//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
354
355void G4PenelopeIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
356 const G4MaterialCutsCouple* couple,
357 const G4DynamicParticle* aDynamicParticle,
358 G4double,
359 G4double)
360{
361 // Penelope model to sample the final state following an hard inelastic interaction.
362 // It makes use of the Generalised Oscillator Strength (GOS) model from
363 // D. Liljequist, J. Phys. D: Appl. Phys. 16 (1983) 1567
364 //
365 // The GOS model is used to calculate the individual cross sections for all
366 // the atomic oscillators coming in the play, taking into account the three
367 // contributions (distant longitudinal collisions, distant transverse collisions and
368 // close collisions). Then the target shell and the interaction channel are
369 // sampled. Final state of the delta-ray (energy, angle) are generated according
370 // to the analytical distributions (dSigma/dW) for the selected interaction
371 // channels.
372 // For e-, the maximum energy for the delta-ray is initialEnergy/2. (because
373 // particles are indistinghusbable), while it is the full initialEnergy for
374 // e+.
375 // The efficiency of the random sampling algorithm (e.g. for close collisions)
376 // decreases when initial and cutoff energy increase (e.g. 87% for 10-keV primary
377 // and 1 keV threshold, 99% for 10-MeV primary and 10-keV threshold).
378 // Differential cross sections have a different form for e+ and e-.
379 //
380 // WARNING: The model provides an _average_ description of inelastic collisions.
381 // Anyway, the energy spectrum associated to distant excitations of a given
382 // atomic shell is approximated as a single resonance. The simulated energy spectra
383 // show _unphysical_ narrow peaks at energies that are multiple of the shell
384 // resonance energies. The spurious speaks are automatically smoothed out after
385 // multiple inelastic collisions.
386 //
387 // The model determines also the original shell from which the delta-ray is expelled,
388 // in order to produce fluorescence de-excitation (from G4DeexcitationManager)
389 //
390 // Fermi density correction is calculated analytically according to
391 // U. Fano, Ann. Rev. Nucl. Sci. 13 (1963),1
392
393 if (verboseLevel > 3)
394 G4cout << "Calling SamplingSecondaries() of G4PenelopeIonisationModel" << G4endl;
395
396 G4double kineticEnergy0 = aDynamicParticle->GetKineticEnergy();
397 const G4ParticleDefinition* theParticle = aDynamicParticle->GetDefinition();
398
399 if (kineticEnergy0 <= LowEnergyLimit())
400 {
401 fParticleChange->SetProposedKineticEnergy(0.);
402 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy0);
403 //Check if there are AtRest processes
404 if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size())
405 //In this case there is at least one AtRest process
406 fParticleChange->ProposeTrackStatus(fStopButAlive);
407 else
408 fParticleChange->ProposeTrackStatus(fStopAndKill);
409 return ;
410 }
411 const G4double electronVolumeDensity =
412 couple->GetMaterial()->GetTotNbOfElectPerVolume();
413 G4ParticleMomentum particleDirection0 = aDynamicParticle->GetMomentumDirection();
414
415 //Initialise final-state variables. The proper values will be set by the methods
416 // CalculateDiscreteForElectrons() and CalculateDiscreteForPositrons()
417 kineticEnergy1=kineticEnergy0;
418 cosThetaPrimary=1.0;
419 energySecondary=0.0;
420 cosThetaSecondary=1.0;
421
422 // Select randomly one element in the current material
423 if (verboseLevel > 2)
424 G4cout << "Going to select element in " << couple->GetMaterial()->GetName() << G4endl;
425
426 //Use sampler of G4CrossSectionHandler for now
427 // atom can be selected effitiantly if element selectors are initialised
428 //const G4Element* anElement = SelectRandomAtom(couple,theParticle,kineticEnergy0);
429
430 G4int iZ = SampleRandomAtom(couple,kineticEnergy0);
431
432 G4double cutForLowEnergySecondaryParticles = 250.0*eV;
433 const G4ProductionCutsTable* theCoupleTable=
434 G4ProductionCutsTable::GetProductionCutsTable();
435 size_t indx = couple->GetIndex();
436 G4double cutE = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
437 G4double cutG = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
438
439 //Production cut for delta-rays (electrons)
440 cutE = std::max(cutForLowEnergySecondaryParticles,cutE);
441 //Production cut for gamma (fluorescence)
442 cutG = std::max(cutForLowEnergySecondaryParticles,cutG);
443
444 if (verboseLevel > 2)
445 G4cout << "Selected Z = " << iZ << G4endl;
446
447 // The method CalculateDiscreteForXXXX() set the private variables:
448 // kineticEnergy1 = energy of the primary electron/positron after the interaction
449 // cosThetaPrimary = std::cos(theta) of the primary after the interaction
450 // energySecondary = energy of the secondary electron
451 // cosThetaSecondary = std::cos(theta) of the secondary
452 if (theParticle == G4Electron::Electron())
453 CalculateDiscreteForElectrons(kineticEnergy0,cutE,iZ,electronVolumeDensity);
454 else if (theParticle == G4Positron::Positron())
455 CalculateDiscreteForPositrons(kineticEnergy0,cutE,iZ,electronVolumeDensity);
456
457 // if (energySecondary == 0) return;
458
459 //Update the primary particle
460 G4double sint = std::sqrt(1. - cosThetaPrimary*cosThetaPrimary);
461 G4double phi = twopi * G4UniformRand();
462 G4double dirx = sint * std::cos(phi);
463 G4double diry = sint * std::sin(phi);
464 G4double dirz = cosThetaPrimary;
465
466 G4ThreeVector electronDirection1(dirx,diry,dirz);
467 electronDirection1.rotateUz(particleDirection0);
468
469 if (kineticEnergy1 > 0)
470 {
471 fParticleChange->ProposeMomentumDirection(electronDirection1);
472 fParticleChange->SetProposedKineticEnergy(kineticEnergy1);
473 }
474 else
475 {
476 fParticleChange->ProposeMomentumDirection(electronDirection1);
477 fParticleChange->SetProposedKineticEnergy(0.*eV);
478 if (theParticle->GetProcessManager()->GetAtRestProcessVector()->size())
479 //In this case there is at least one AtRest process
480 fParticleChange->ProposeTrackStatus(fStopButAlive);
481 else
482 fParticleChange->ProposeTrackStatus(fStopAndKill);
483 }
484
485 //Generate the delta ray
486 G4int iosc2 = 0;
487 G4double ioniEnergy = 0.0;
488 if (iOsc > 0) {
489 ioniEnergy=(*(ionizationEnergy->find(iZ)->second))[iOsc];
490 iosc2 = (ionizationEnergy->find(iZ)->second->size()) - iOsc; //they are in reversed order
491 }
492
493 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
494 G4double bindingEnergy = 0.0;
495 G4int shellId = 0;
496 if (iOsc > 0)
497 {
498 const G4AtomicShell* shell = transitionManager->Shell(iZ,iosc2-1); // Modified by Alf
499 bindingEnergy = shell->BindingEnergy();
500 shellId = shell->ShellId();
501 }
502
503 G4double ionEnergy = bindingEnergy; //energy spent to ionise the atom according to G4dabatase
504 G4double eKineticEnergy = energySecondary;
505
506 //This is an awful thing: Penelope generates the fluorescence only for L and K shells
507 //(i.e. Osc = 1 --> 4). For high-Z, the other shells can be quite relevant. In this case
508 //one MUST ensure ''by hand'' the energy conservation. Then there is the other problem that
509 //the fluorescence database of Penelope doesn not match that of Geant4.
510
511 G4double energyBalance = kineticEnergy0 - kineticEnergy1 - energySecondary; //Penelope Balance
512
513 if (std::abs(energyBalance) < 1*eV)
514 //in this case Penelope didn't subtract the fluorescence energy: do here by hand
515 eKineticEnergy = energySecondary - bindingEnergy;
516 else
517 //Penelope subtracted the fluorescence, but one has to match the databases
518 eKineticEnergy = energySecondary+ioniEnergy-bindingEnergy;
[1007]519
520 //VERIFICA QUI LA STORIA DEL LOCAL ENERGY DEPOSIT!
521
[968]522 G4double localEnergyDeposit = ionEnergy;
523 G4double energyInFluorescence = 0.0*eV;
524
525 std::vector<G4DynamicParticle*> *photonVector = 0;
526 if (fUseAtomicDeexcitation)
527 {
528 if (iZ>5 && (ionEnergy > cutG || ionEnergy > cutE))
529 {
530 photonVector = deexcitationManager.GenerateParticles(iZ,shellId);
531 //Check for single photons if they are above threshold
532 for (size_t k=0;k<photonVector->size();k++)
533 {
534 G4DynamicParticle* aPhoton = (*photonVector)[k];
535 if (aPhoton)
536 {
537 G4double itsCut = cutG;
538 if (aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
539 G4double itsEnergy = aPhoton->GetKineticEnergy();
540 if (itsEnergy > itsCut && itsEnergy <= ionEnergy)
541 {
542 localEnergyDeposit -= itsEnergy;
543 energyInFluorescence += itsEnergy;
544 }
545 else
546 {
547 delete aPhoton;
548 (*photonVector)[k] = 0;
549 }
550 }
551 }
552 }
553 }
554
555 // Generate the delta ray
556 G4double sin2 = std::sqrt(1. - cosThetaSecondary*cosThetaSecondary);
557 G4double phi2 = twopi * G4UniformRand();
558
559 G4double xEl = sin2 * std::cos(phi2);
560 G4double yEl = sin2 * std::sin(phi2);
561 G4double zEl = cosThetaSecondary;
562 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
563 eDirection.rotateUz(particleDirection0);
564
565 G4DynamicParticle* deltaElectron = new G4DynamicParticle (G4Electron::Electron(),
566 eDirection,eKineticEnergy) ;
567 fvect->push_back(deltaElectron);
568
569 //Generate fluorescence, if it is the case
570 //This block is executed only if there is at least one secondary photon produced by
571 //G4AtomicDeexcitation
572 if (photonVector)
573 {
574 for (size_t ll=0;ll<photonVector->size();ll++)
575 if ((*photonVector)[ll])
576 {
577 G4DynamicParticle* aFluorescencePhoton = (*photonVector)[ll];
578 fvect->push_back(aFluorescencePhoton);
579 }
580 }
581 delete photonVector;
582
583 if (localEnergyDeposit < 0)
584 {
585 G4cout << "WARNING-"
586 << "G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
587 << G4endl;
588 localEnergyDeposit=0.;
589 }
590 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
591
592 if (verboseLevel > 1)
593 {
594 G4cout << "-----------------------------------------------------------" << G4endl;
595 G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
596 G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
597 G4cout << "-----------------------------------------------------------" << G4endl;
598 G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
599 G4cout << "Delta ray " << eKineticEnergy/keV << " keV" << G4endl;
600 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
601 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
602 G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+kineticEnergy1+
603 localEnergyDeposit)/keV <<
604 " keV" << G4endl;
605 G4cout << "-----------------------------------------------------------" << G4endl;
606 }
607 if (verboseLevel > 0)
608 {
609 G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+kineticEnergy1+
610 localEnergyDeposit-kineticEnergy0);
611 if (energyDiff > 0.05*keV)
612 G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<
613 (eKineticEnergy+energyInFluorescence+kineticEnergy1+localEnergyDeposit)/keV <<
614 " keV (final) vs. " <<
615 kineticEnergy0/keV << " keV (initial)" << G4endl;
616 }
617}
618
619//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
620
621void G4PenelopeIonisationModel::ReadData()
622{
623 if (verboseLevel > 2)
624 {
625 G4cout << "Data from G4PenelopeIonisationModel read " << G4endl;
626 }
627 char* path = getenv("G4LEDATA");
628 if (!path)
629 {
630 G4String excep = "G4PenelopeIonisationModel - G4LEDATA environment variable not set!";
631 G4Exception(excep);
632 }
633 G4String pathString(path);
634 G4String pathFile = pathString + "/penelope/ion-pen.dat";
635 std::ifstream file(pathFile);
636
637 if (!file.is_open())
638 {
639 G4String excep = "G4PenelopeIonisationModel - data file " + pathFile + " not found!";
640 G4Exception(excep);
641 }
642
643 if (!ionizationEnergy || !resonanceEnergy || !occupationNumber || !shellFlag)
644 {
645 G4String excep = "G4PenelopeIonisationModel: problem with reading data from file";
646 G4Exception(excep);
647 }
648
649 G4int Z=1,nLevels=0;
650 G4int test,test1;
651
652 do{
653 G4DataVector* occVector = new G4DataVector;
654 G4DataVector* ionEVector = new G4DataVector;
655 G4DataVector* resEVector = new G4DataVector;
656 G4DataVector* shellIndVector = new G4DataVector;
657 file >> Z >> nLevels;
658 G4double a1,a2,a3,a4;
659 G4int k1,k2,k3;
660 for (G4int h=0;h<nLevels;h++)
661 {
662 //index,occup number,ion energy,res energy,fj0,kz,shell flag
663 file >> k1 >> a1 >> a2 >> a3 >> a4 >> k2 >> k3;
664 //Make explicit unit of measurements for ionisation and resonance energy,
665 // which is MeV
666 a2 *= MeV;
667 a3 *= MeV;
668 //
669 occVector->push_back(a1);
670 ionEVector->push_back(a2);
671 resEVector->push_back(a3);
672 shellIndVector->push_back((G4double) k3);
673 }
674 //Ok, done for element Z
675 occupationNumber->insert(std::make_pair(Z,occVector));
676 ionizationEnergy->insert(std::make_pair(Z,ionEVector));
677 resonanceEnergy->insert(std::make_pair(Z,resEVector));
678 shellFlag->insert(std::make_pair(Z,shellIndVector));
679 file >> test >> test1; //-1 -1 close the data for each Z
680 if (test > 0) {
681 G4String excep = "G4PenelopeIonisationModel - data file corrupted!";
682 G4Exception(excep);
683 }
684 }while (test != -2); //the very last Z is closed with -2 instead of -1
685}
686
687
688//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
689
690G4double G4PenelopeIonisationModel::CalculateDeltaFermi(G4double kinEnergy ,G4int Z,
691 G4double electronVolumeDensity)
692{
693 G4double plasmaEnergyCoefficient = 1.377e-39*(MeV*MeV*m3); //(e*hbar)^2/(epsilon0*electron_mass)
694 G4double plasmaEnergySquared = plasmaEnergyCoefficient*electronVolumeDensity;
695 // std::sqrt(plasmaEnergySquared) is the plasma energy of the solid (MeV)
696 G4double gam = 1.0+kinEnergy/electron_mass_c2;
697 G4double gam2=gam*gam;
698 G4double delta = 0.0;
699
700 //Density effect
701 G4double TST = ((G4double) Z)/(gam2*plasmaEnergySquared);
702
703 G4double wl2 = 0.0;
704 G4double fdel=0.0;
705 G4double wr=0;
706 size_t nbOsc = resonanceEnergy->find(Z)->second->size();
707 for(size_t i=0;i<nbOsc;i++)
708 {
709 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
710 wr = (*(resonanceEnergy->find(Z)->second))[i];
711 fdel += occupNb/(wr*wr+wl2);
712 }
713 if (fdel < TST) return delta;
714 G4double help1 = (*(resonanceEnergy->find(Z)->second))[nbOsc-1];
715 wl2 = help1*help1;
716 do{
717 wl2=wl2*2.0;
718 fdel = 0.0;
719 for (size_t ii=0;ii<nbOsc;ii++){
720 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[ii];
721 wr = (*(resonanceEnergy->find(Z)->second))[ii];
722 fdel += occupNb/(wr*wr+wl2);
723 }
724 }while (fdel > TST);
725 G4double wl2l=0.0;
726 G4double wl2u = wl2;
727 G4double control = 0.0;
728 do{
729 wl2=0.5*(wl2l+wl2u);
730 fdel = 0.0;
731 for (size_t jj=0;jj<nbOsc;jj++){
732 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[jj];
733 wr = (*(resonanceEnergy->find(Z)->second))[jj];
734 fdel += occupNb/(wr*wr+wl2);
735 }
736 if (fdel > TST)
737 wl2l = wl2;
738 else
739 wl2u = wl2;
740 control = wl2u-wl2l-wl2*1e-12;
741 }while(control>0);
742
743 //Density correction effect
744 for (size_t kk=0;kk<nbOsc;kk++){
745 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[kk];
746 wr = (*(resonanceEnergy->find(Z)->second))[kk];
747 delta += occupNb*std::log(1.0+wl2/(wr*wr));
748 }
749 delta = (delta/((G4double) Z))-wl2/(gam2*plasmaEnergySquared);
750 return delta;
751}
752
753//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
754
755void G4PenelopeIonisationModel::CalculateDiscreteForElectrons(G4double kinEnergy,G4double cutoffEnergy,
756 G4int Z,G4double electronVolumeDensity)
757{
758 if (verboseLevel > 2)
759 G4cout << "Entering in CalculateDiscreteForElectrons() for energy " <<
760 kinEnergy/keV << " keV " << G4endl;
761
762 //Initialization of variables to be calculated in this routine
763 kineticEnergy1=kinEnergy;
764 cosThetaPrimary=1.0;
765 energySecondary=0.0;
766 cosThetaSecondary=1.0;
767 iOsc=-1;
768
769 //constants
770 G4double rb=kinEnergy+2.0*electron_mass_c2;
771 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
772 G4double gamma2 = gamma*gamma;
773 G4double beta2 = (gamma2-1.0)/gamma2;
774 G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
775 G4double cps = kinEnergy*rb;
776 G4double cp = std::sqrt(cps);
777
778 G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
779 G4double distantTransvCS0 = std::max(std::log(gamma2)-beta2-delta,0.0);
780
781 G4double rl,rl1;
782
783 if (cutoffEnergy > kinEnergy) return; //delta rays are not generated
784
785 G4DataVector* qm = new G4DataVector();
786 G4DataVector* cumulHardCS = new G4DataVector();
787 std::vector<G4int> *typeOfInteraction = new std::vector<G4int>;
788 G4DataVector* nbOfLevel = new G4DataVector();
789
790 //Hard close collisions with outer shells
791 G4double wmaxc = 0.5*kinEnergy;
792 G4double closeCS0 = 0.0;
793 G4double closeCS = 0.0;
794 if (cutoffEnergy>0.1*eV)
795 {
796 rl=cutoffEnergy/kinEnergy;
797 rl1=1.0-rl;
798 if (rl < 0.5)
799 closeCS0 = (amol*(0.5-rl)+(1.0/rl)-(1.0/rl1)+(1.0-amol)*std::log(rl/rl1))/kinEnergy;
800 }
801
802 // Cross sections for the different oscillators
803
804 // totalHardCS contains the cumulative hard interaction cross section for the different
805 // excitable levels and the different interaction channels (close, distant, etc.),
806 // i.e.
807 // cumulHardCS[0] = 0.0
808 // cumulHardCS[1] = 1st excitable level (distant longitudinal only)
809 // cumulHardCS[2] = 1st excitable level (distant longitudinal + transverse)
810 // cumulHardCS[3] = 1st excitable level (distant longitudinal + transverse + close)
811 // cumulHardCS[4] = 1st excitable level (all channels) + 2nd excitable level (distant long only)
812 // etc.
813 // This is used for sampling the atomic level which is ionised and the channel of the
814 // interaction.
815 //
816 // For each index iFill of the cumulHardCS vector,
817 // nbOfLevel[iFill] contains the current excitable atomic level and
818 // typeOfInteraction[iFill] contains the current interaction channel, with the legenda:
819 // 1 = distant longitudinal interaction
820 // 2 = distant transverse interaction
821 // 3 = close collision
822 // 4 = close collision with outer shells (in this case nbOfLevel < 0 --> no binding energy)
823
824
825 G4int nOscil = ionizationEnergy->find(Z)->second->size();
826 G4double totalHardCS = 0.0;
827 G4double involvedElectrons = 0.0;
828 for (G4int i=0;i<nOscil;i++){
829 G4double wi = (*(resonanceEnergy->find(Z)->second))[i];
830 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
831
832 //Distant excitations
833 if (wi>cutoffEnergy && wi<kinEnergy)
834 {
835 if (wi>(1e-6*kinEnergy))
836 {
837 G4double cpp=std::sqrt((kinEnergy-wi)*(kinEnergy-wi+2.0*electron_mass_c2));
838 qm->push_back(std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2);
839 }
840 else
841 {
842 qm->push_back((wi*wi)/(beta2*2.0*electron_mass_c2));
843 }
844 if ((*qm)[i] < wi)
845 {
846
847 G4double distantLongitCS = occupNb*std::log(wi*((*qm)[i]+2.0*electron_mass_c2)/
848 ((*qm)[i]*(wi+2.0*electron_mass_c2)))/wi;
849 cumulHardCS->push_back(totalHardCS);
850 typeOfInteraction->push_back(1); //distant longitudinal
851 nbOfLevel->push_back((G4double) i); //only excitable level are counted
852 totalHardCS += distantLongitCS;
853
854 G4double distantTransvCS = occupNb*distantTransvCS0/wi;
855
856 cumulHardCS->push_back(totalHardCS);
857 typeOfInteraction->push_back(2); //distant tranverse
858 nbOfLevel->push_back((G4double) i);
859 totalHardCS += distantTransvCS;
860 }
861 }
862 else
863 qm->push_back(wi);
864
865 //close collisions
866 if(wi < wmaxc)
867 {
868 if (wi < cutoffEnergy)
869 involvedElectrons += occupNb;
870 else
871 {
872 rl=wi/kinEnergy;
873 rl1=1.0-rl;
874 closeCS = occupNb*(amol*(0.5-rl)+(1.0/rl)-(1.0/rl1)+(1.0-amol)*std::log(rl/rl1))/kinEnergy;
875 cumulHardCS->push_back(totalHardCS);
876 typeOfInteraction->push_back(3); //close
877 nbOfLevel->push_back((G4double) i);
878 totalHardCS += closeCS;
879 }
880 }
881 } // loop on the levels
882
883 cumulHardCS->push_back(totalHardCS);
884 typeOfInteraction->push_back(4); //close interaction with outer shells
885 nbOfLevel->push_back(-1.0);
886 totalHardCS += involvedElectrons*closeCS0;
887 cumulHardCS->push_back(totalHardCS); //this is the final value of the totalHardCS
888
889 if (totalHardCS < 1e-30) {
890 kineticEnergy1=kinEnergy;
891 cosThetaPrimary=1.0;
892 energySecondary=0.0;
893 cosThetaSecondary=0.0;
894 iOsc=-1;
895 delete qm;
896 delete cumulHardCS;
897 delete typeOfInteraction;
898 delete nbOfLevel;
899 return;
900 }
901
902 //Testing purposes
903 if (verboseLevel > 6)
904 {
905 for (size_t t=0;t<cumulHardCS->size();t++)
906 G4cout << (*cumulHardCS)[t] << " " << (*typeOfInteraction)[t] <<
907 " " << (*nbOfLevel)[t] << G4endl;
908 }
909
910 //Selection of the active oscillator on the basis of the cumulative cross sections
911 G4double TST = totalHardCS*G4UniformRand();
912 G4int is=0;
913 G4int js= nbOfLevel->size();
914 do{
915 G4int it=(is+js)/2;
916 if (TST > (*cumulHardCS)[it]) is=it;
917 if (TST <= (*cumulHardCS)[it]) js=it;
918 }while((js-is) > 1);
919
920 G4double UII=0.0;
921 G4double rkc=cutoffEnergy/kinEnergy;
922 G4double dde;
923 G4int kks;
924
925 G4int sampledInteraction = (*typeOfInteraction)[is];
926 iOsc = (G4int) (*nbOfLevel)[is];
927
928 if (verboseLevel > 4)
929 G4cout << "Chosen interaction #:" << sampledInteraction << " on level " << iOsc << G4endl;
930
931 //Generates the final state according to the sampled level and
932 //interaction channel
933
934 if (sampledInteraction == 1) //Hard distant longitudinal collisions
935 {
936 dde= (*(resonanceEnergy->find(Z)->second))[iOsc];
937 kineticEnergy1=kinEnergy-dde;
938 G4double qs=(*qm)[iOsc]/(1.0+((*qm)[iOsc]/(2.0*electron_mass_c2)));
939 G4double q=qs/(std::pow((qs/dde)*(1.0+(0.5*dde/electron_mass_c2)),G4UniformRand())-(0.5*qs/electron_mass_c2));
940 G4double qtrev = q*(q+2.0*electron_mass_c2);
941 G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
942 cosThetaPrimary = (cpps+cps-qtrev)/(2.0*cp*std::sqrt(cpps));
943 if (cosThetaPrimary>1.0) cosThetaPrimary=1.0;
944 //Energy and emission angle of the delta ray
945 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
946 //kks > 4 means that we are in an outer shell
947 if (kks>4)
948 energySecondary=dde;
949 else
950 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
951
952 cosThetaSecondary = 0.5*(dde*(kinEnergy+rb-dde)+qtrev)/std::sqrt(cps*qtrev);
953 if (cosThetaSecondary>1.0) cosThetaSecondary=1.0;
954 }
955 else if (sampledInteraction == 2) //Hard distant transverse collisions
956 {
957 dde=(*(resonanceEnergy->find(Z)->second))[iOsc];
958 kineticEnergy1=kinEnergy-dde;
959 cosThetaPrimary=1.0;
960 //Energy and emission angle of the delta ray
961 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
962 if (kks>4)
963 {
964 energySecondary=dde;
965 }
966 else
967 {
968 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
969 }
970 cosThetaSecondary = 1.0;
971 }
972 else if (sampledInteraction == 3 || sampledInteraction == 4) //Close interaction
973 {
974 if (sampledInteraction == 4) //interaction with inner shells
975 {
976 UII=0.0;
977 rkc = cutoffEnergy/kinEnergy;
978 iOsc = -1;
979 }
980 else
981 {
982 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
983 if (kks > 4) {
984 UII=0.0;
985 }
986 else
987 {
988 UII = (*(ionizationEnergy->find(Z)->second))[iOsc];
989 }
990 rkc = (*(resonanceEnergy->find(Z)->second))[iOsc]/kinEnergy;
991 }
992 G4double A = 0.5*amol;
993 G4double arkc = A*0.5*rkc;
994 G4double phi,rk2,rk,rkf;
995 do{
996 G4double fb = (1.0+arkc)*G4UniformRand();
997 if (fb<1.0)
998 {
999 rk=rkc/(1.0-fb*(1.0-(rkc*2.0)));
1000 }
1001 else{
1002 rk = rkc+(fb-1.0)*(0.5-rkc)/arkc;
1003 }
1004 rk2 = rk*rk;
1005 rkf = rk/(1.0-rk);
1006 phi = 1.0+(rkf*rkf)-rkf+amol*(rk2+rkf);
1007 }while ((G4UniformRand()*(1.0+A*rk2)) > phi);
1008 //Energy and scattering angle (primary electron);
1009 kineticEnergy1 = kinEnergy*(1.0-rk);
1010 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kinEnergy*(rb-(rk*kinEnergy))));
1011 //Energy and scattering angle of the delta ray
1012 energySecondary = kinEnergy-kineticEnergy1-UII;
1013 cosThetaSecondary = std::sqrt(rk*kinEnergy*rb/(kinEnergy*(rk*kinEnergy+2.0*electron_mass_c2)));
1014 }
1015 else
1016 {
1017 G4String excep = "G4PenelopeIonisationModel - Error in the calculation of the final state";
1018 G4Exception(excep);
1019 }
1020
1021 delete qm;
1022 delete cumulHardCS;
1023
1024 delete typeOfInteraction;
1025 delete nbOfLevel;
1026
1027 return;
1028}
1029
1030//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1031
1032 void G4PenelopeIonisationModel::CalculateDiscreteForPositrons(G4double kinEnergy,G4double cutoffEnergy,
1033 G4int Z,G4double electronVolumeDensity)
1034{
1035 kineticEnergy1=kinEnergy;
1036 cosThetaPrimary=1.0;
1037 energySecondary=0.0;
1038 cosThetaSecondary=1.0;
1039 iOsc=-1;
1040 //constants
1041 G4double rb=kinEnergy+2.0*electron_mass_c2;
1042 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1043 G4double gamma2 = gamma*gamma;
1044 G4double beta2 = (gamma2-1.0)/gamma2;
1045 G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
1046 G4double cps = kinEnergy*rb;
1047 G4double cp = std::sqrt(cps);
1048 G4double help = (gamma+1.0)*(gamma+1.0);
1049 G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1050 G4double bha2 = amol*(3.0+1.0/help);
1051 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1052 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1053
1054 G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
1055 G4double distantTransvCS0 = std::max(std::log(gamma2)-beta2-delta,0.0);
1056
1057 G4double rl,rl1;
1058
1059 if (cutoffEnergy > kinEnergy) return; //delta rays are not generated
1060
1061 G4DataVector* qm = new G4DataVector();
1062 G4DataVector* cumulHardCS = new G4DataVector();
1063 std::vector<G4int> *typeOfInteraction = new std::vector<G4int>;
1064 G4DataVector* nbOfLevel = new G4DataVector();
1065
1066
1067 //Hard close collisions with outer shells
1068 G4double wmaxc = kinEnergy;
1069 G4double closeCS0 = 0.0;
1070 G4double closeCS = 0.0;
1071 if (cutoffEnergy>0.1*eV)
1072 {
1073 rl=cutoffEnergy/kinEnergy;
1074 rl1=1.0-rl;
1075 if (rl < 1.0)
1076 closeCS0 = (((1.0/rl)-1.0) + bha1*std::log(rl) + bha2*rl1
1077 + (bha3/2.0)*((rl*rl)-1.0)
1078 + (bha4/3.0)*(1.0-(rl*rl*rl)))/kinEnergy;
1079 }
1080
1081 // Cross sections for the different oscillators
1082
1083 // totalHardCS contains the cumulative hard interaction cross section for the different
1084 // excitable levels and the different interaction channels (close, distant, etc.),
1085 // i.e.
1086 // cumulHardCS[0] = 0.0
1087 // cumulHardCS[1] = 1st excitable level (distant longitudinal only)
1088 // cumulHardCS[2] = 1st excitable level (distant longitudinal + transverse)
1089 // cumulHardCS[3] = 1st excitable level (distant longitudinal + transverse + close)
1090 // cumulHardCS[4] = 1st excitable level (all channels) + 2nd excitable level (distant long only)
1091 // etc.
1092 // This is used for sampling the atomic level which is ionised and the channel of the
1093 // interaction.
1094 //
1095 // For each index iFill of the cumulHardCS vector,
1096 // nbOfLevel[iFill] contains the current excitable atomic level and
1097 // typeOfInteraction[iFill] contains the current interaction channel, with the legenda:
1098 // 1 = distant longitudinal interaction
1099 // 2 = distant transverse interaction
1100 // 3 = close collision
1101 // 4 = close collision with outer shells (in this case nbOfLevel < 0 --> no binding energy)
1102
1103
1104 G4int nOscil = ionizationEnergy->find(Z)->second->size();
1105 G4double totalHardCS = 0.0;
1106 G4double involvedElectrons = 0.0;
1107 for (G4int i=0;i<nOscil;i++){
1108 G4double wi = (*(resonanceEnergy->find(Z)->second))[i];
1109 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
1110 //Distant excitations
1111 if (wi>cutoffEnergy && wi<kinEnergy)
1112 {
1113 if (wi>(1e-6*kinEnergy)){
1114 G4double cpp=std::sqrt((kinEnergy-wi)*(kinEnergy-wi+2.0*electron_mass_c2));
1115 qm->push_back(std::sqrt((cp-cpp)*(cp-cpp)+ electron_mass_c2 * electron_mass_c2)-electron_mass_c2);
1116 }
1117 else
1118 {
1119 qm->push_back(wi*wi/(beta2+2.0*electron_mass_c2));
1120 }
1121 //verificare che quando arriva qui il vettore ha SEMPRE l'i-esimo elemento
1122 if ((*qm)[i] < wi)
1123 {
1124
1125 G4double distantLongitCS = occupNb*std::log(wi*((*qm)[i]+2.0*electron_mass_c2)/
1126 ((*qm)[i]*(wi+2.0*electron_mass_c2)))/wi;
1127 cumulHardCS->push_back(totalHardCS);
1128 typeOfInteraction->push_back(1); //distant longitudinal
1129 nbOfLevel->push_back((G4double) i); //only excitable level are counted
1130 totalHardCS += distantLongitCS;
1131
1132 G4double distantTransvCS = occupNb*distantTransvCS0/wi;
1133
1134 cumulHardCS->push_back(totalHardCS);
1135 typeOfInteraction->push_back(2); //distant tranverse
1136 nbOfLevel->push_back((G4double) i);
1137 totalHardCS += distantTransvCS;
1138 }
1139 }
1140 else
1141 qm->push_back(wi);
1142
1143 //close collisions
1144 if(wi < wmaxc)
1145 {
1146 if (wi < cutoffEnergy) {
1147 involvedElectrons += occupNb;
1148 }
1149 else
1150 {
1151 rl=wi/kinEnergy;
1152 rl1=1.0-rl;
1153 closeCS = occupNb*(((1.0/rl)-1.0)+bha1*std::log(rl)+bha2*rl1
1154 + (bha3/2.0)*((rl*rl)-1.0)
1155 + (bha4/3.0)*(1.0-(rl*rl*rl)))/kinEnergy;
1156 cumulHardCS->push_back(totalHardCS);
1157 typeOfInteraction->push_back(3); //close
1158 nbOfLevel->push_back((G4double) i);
1159 totalHardCS += closeCS;
1160 }
1161 }
1162 } // loop on the levels
1163
1164 cumulHardCS->push_back(totalHardCS);
1165 typeOfInteraction->push_back(4); //close interaction with outer shells
1166 nbOfLevel->push_back(-1.0);
1167 totalHardCS += involvedElectrons*closeCS0;
1168 cumulHardCS->push_back(totalHardCS); //this is the final value of the totalHardCS
1169
1170 if (totalHardCS < 1e-30) {
1171 kineticEnergy1=kinEnergy;
1172 cosThetaPrimary=1.0;
1173 energySecondary=0.0;
1174 cosThetaSecondary=0.0;
1175 iOsc=-1;
1176 delete qm;
1177 delete cumulHardCS;
1178 delete typeOfInteraction;
1179 delete nbOfLevel;
1180 return;
1181 }
1182
1183 //Selection of the active oscillator on the basis of the cumulative cross sections
1184 G4double TST = totalHardCS*G4UniformRand();
1185 G4int is=0;
1186 G4int js= nbOfLevel->size();
1187 do{
1188 G4int it=(is+js)/2;
1189 if (TST > (*cumulHardCS)[it]) is=it;
1190 if (TST <= (*cumulHardCS)[it]) js=it;
1191 }while((js-is) > 1);
1192
1193 G4double UII=0.0;
1194 G4double rkc=cutoffEnergy/kinEnergy;
1195 G4double dde;
1196 G4int kks;
1197
1198 G4int sampledInteraction = (*typeOfInteraction)[is];
1199 iOsc = (G4int) (*nbOfLevel)[is];
1200
1201 //Generates the final state according to the sampled level and
1202 //interaction channel
1203
1204 if (sampledInteraction == 1) //Hard distant longitudinal collisions
1205 {
1206 dde= (*(resonanceEnergy->find(Z)->second))[iOsc];
1207 kineticEnergy1=kinEnergy-dde;
1208 G4double qs=(*qm)[iOsc]/(1.0+((*qm)[iOsc]/(2.0*electron_mass_c2)));
1209 G4double q=qs/(std::pow((qs/dde)*(1.0+(0.5*dde/electron_mass_c2)),G4UniformRand())-(0.5*qs/electron_mass_c2));
1210 G4double qtrev = q*(q+2.0*electron_mass_c2);
1211 G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
1212 cosThetaPrimary = (cpps+cps-qtrev)/(2.0*cp*std::sqrt(cpps));
1213 if (cosThetaPrimary>1.0) cosThetaPrimary=1.0;
1214 //Energy and emission angle of the delta ray
1215 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1216 if (kks>4)
1217 energySecondary=dde;
1218
1219 else
1220 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
1221 cosThetaSecondary = 0.5*(dde*(kinEnergy+rb-dde)+qtrev)/std::sqrt(cps*qtrev);
1222 if (cosThetaSecondary>1.0) cosThetaSecondary=1.0;
1223 }
1224 else if (sampledInteraction == 2) //Hard distant transverse collisions
1225 {
1226 dde=(*(resonanceEnergy->find(Z)->second))[iOsc];
1227 kineticEnergy1=kinEnergy-dde;
1228 cosThetaPrimary=1.0;
1229 //Energy and emission angle of the delta ray
1230 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1231 if (kks>4)
1232 {
1233 energySecondary=dde;
1234 }
1235 else
1236 {
1237 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
1238 }
1239 cosThetaSecondary = 1.0;
1240 }
1241 else if (sampledInteraction == 3 || sampledInteraction == 4) //Close interaction
1242 {
1243 if (sampledInteraction == 4) //interaction with inner shells
1244 {
1245 UII=0.0;
1246 rkc = cutoffEnergy/kinEnergy;
1247 iOsc = -1;
1248 }
1249 else
1250 {
1251 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1252 if (kks > 4) {
1253 UII=0.0;
1254 }
1255 else
1256 {
1257 UII = (*(ionizationEnergy->find(Z)->second))[iOsc];
1258 }
1259 rkc = (*(resonanceEnergy->find(Z)->second))[iOsc]/kinEnergy;
1260 }
1261 G4double phi,rk;
1262 do{
1263 rk=rkc/(1.0-G4UniformRand()*(1.0-rkc));
1264 phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1265 }while ( G4UniformRand() > phi);
1266 //Energy and scattering angle (primary electron);
1267 kineticEnergy1 = kinEnergy*(1.0-rk);
1268 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kinEnergy*(rb-(rk*kinEnergy))));
1269 //Energy and scattering angle of the delta ray
1270 energySecondary = kinEnergy-kineticEnergy1-UII;
1271 cosThetaSecondary = std::sqrt(rk*kinEnergy*rb/(kinEnergy*(rk*kinEnergy+2.0*electron_mass_c2)));
1272 }
1273 else
1274 {
1275 G4String excep = "G4PenelopeIonisationModel - Error in the calculation of the final state";
1276 G4Exception(excep);
1277 }
1278
1279 delete qm;
1280 delete cumulHardCS;
1281 delete typeOfInteraction;
1282 delete nbOfLevel;
1283
1284 return;
1285}
1286
1287//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1288
1289G4double G4PenelopeIonisationModel::CalculateCrossSectionsRatio(G4double kinEnergy,
1290 G4double cutoffEnergy,
1291 G4int Z, G4double electronVolumeDensity,
1292 const G4ParticleDefinition* theParticle)
1293{
1294 //Constants
1295 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1296 G4double gamma2 = gamma*gamma;
1297 G4double beta2 = (gamma2-1.0)/gamma2;
1298 G4double constant = pi*classic_electr_radius*classic_electr_radius*2.0*electron_mass_c2/beta2;
1299 G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
1300 G4int nbOsc = (G4int) resonanceEnergy->find(Z)->second->size();
1301 G4double softCS = 0.0;
1302 G4double hardCS = 0.0;
1303 for (G4int i=0;i<nbOsc;i++){
1304 G4double resEnergy = (*(resonanceEnergy->find(Z)->second))[i];
1305 std::pair<G4double,G4double> SoftAndHardXS(0.,0.);
1306 if (theParticle == G4Electron::Electron())
1307 SoftAndHardXS = CrossSectionsRatioForElectrons(kinEnergy,resEnergy,delta,cutoffEnergy);
1308 else if (theParticle == G4Positron::Positron())
1309 SoftAndHardXS = CrossSectionsRatioForPositrons(kinEnergy,resEnergy,delta,cutoffEnergy);
1310 G4double occupNb = (*(occupationNumber->find(Z)->second))[i];
1311 softCS += occupNb*constant*SoftAndHardXS.first;
1312 hardCS += occupNb*constant*SoftAndHardXS.second;
1313 }
1314 G4double ratio = 0.0;
1315
1316 if (softCS+hardCS) ratio = (hardCS)/(softCS+hardCS);
1317 return ratio;
1318}
1319
1320//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1321
1322std::pair<G4double,G4double> G4PenelopeIonisationModel::CrossSectionsRatioForElectrons(G4double kineticEnergy,
1323 G4double resEnergy,
1324 G4double densityCorrection,
1325 G4double cutoffEnergy)
1326{
1327 std::pair<G4double,G4double> theResult(0.,0.);
1328 if (kineticEnergy < resEnergy) return theResult;
1329
1330 //Calculate constants
1331 G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
1332 G4double gamma2 = gamma*gamma;
1333 G4double beta2 = (gamma2-1.0)/gamma2;
1334 G4double cps = kineticEnergy*(kineticEnergy+2.0*electron_mass_c2);
1335 G4double amol = (kineticEnergy/(kineticEnergy+electron_mass_c2)) * (kineticEnergy/(kineticEnergy+electron_mass_c2)) ;
1336 G4double hardCont = 0.0;
1337 G4double softCont = 0.0;
1338
1339 //Distant interactions
1340 G4double cp1s = (kineticEnergy-resEnergy)*(kineticEnergy-resEnergy+2.0*electron_mass_c2);
1341 G4double cp1 = std::sqrt(cp1s);
1342 G4double cp = std::sqrt(cps);
1343 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1344
1345 //Distant longitudinal interactions
1346 G4double qm = 0.0;
1347
1348 if (resEnergy > kineticEnergy*(1e-6))
1349 {
1350 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1351 }
1352 else
1353 {
1354 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1355 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1356 }
1357
1358 if (qm < resEnergy)
1359 {
1360 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1361 }
1362 else
1363 {
1364 sdLong = 0.0;
1365 }
1366
1367 if (sdLong > 0) {
1368 sdTrans = std::max(std::log(gamma2)-beta2-densityCorrection,0.0);
1369 sdDist = sdTrans + sdLong;
1370 if (cutoffEnergy > resEnergy)
1371 {
1372 softCont = sdDist/resEnergy;
1373 }
1374 else
1375 {
1376 hardCont = sdDist/resEnergy;
1377 }
1378 }
1379
1380 // Close collisions (Moeller's cross section)
1381 G4double wl = std::max(cutoffEnergy,resEnergy);
1382 G4double wu = 0.5*kineticEnergy;
1383
1384 if (wl < (wu-1*eV))
1385 {
1386 hardCont += (1.0/(kineticEnergy-wu))-(1.0/(kineticEnergy-wl))
1387 - (1.0/wu)+(1.0/wl)
1388 + (1.0-amol)*std::log(((kineticEnergy-wu)*wl)/((kineticEnergy-wl)*wu))/kineticEnergy
1389 + amol*(wu-wl)/(kineticEnergy*kineticEnergy);
1390 wu=wl;
1391 }
1392
1393 wl = resEnergy;
1394 if (wl > (wu-1*eV))
1395 {
1396 theResult.first = softCont;
1397 theResult.second = hardCont;
1398 return theResult;
1399 }
1400 softCont += (1.0/(kineticEnergy-wu))-(1.0/(kineticEnergy-wl))
1401 - (1.0/wu)+(1.0/wl)
1402 + (1.0-amol)*std::log(((kineticEnergy-wu)*wl)/((kineticEnergy-wl)*wu))/kineticEnergy
1403 + amol*(wu-wl)/(kineticEnergy*kineticEnergy);
1404 theResult.first = softCont;
1405 theResult.second = hardCont;
1406 return theResult;
1407}
1408
1409
1410//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1411
1412std::pair<G4double,G4double> G4PenelopeIonisationModel::CrossSectionsRatioForPositrons(G4double kineticEnergy,
1413 G4double resEnergy,
1414 G4double densityCorrection,
1415 G4double cutoffEnergy)
1416{
1417
1418 std::pair<G4double,G4double> theResult(0.,0.);
1419 if (kineticEnergy < resEnergy) return theResult;
1420
1421 //Calculate constants
1422 G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
1423 G4double gamma2 = gamma*gamma;
1424 G4double beta2 = (gamma2-1.0)/gamma2;
1425 G4double cps = kineticEnergy*(kineticEnergy+2.0*electron_mass_c2);
1426 G4double amol = (kineticEnergy/(kineticEnergy+electron_mass_c2)) * (kineticEnergy/(kineticEnergy+electron_mass_c2)) ;
1427 G4double help = (gamma+1.0)*(gamma+1.0);
1428 G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1429 G4double bha2 = amol*(3.0+1.0/help);
1430 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1431 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1432 G4double hardCont = 0.0;
1433 G4double softCont = 0.0;
1434
1435 //Distant interactions
1436 G4double cp1s = (kineticEnergy-resEnergy)*(kineticEnergy-resEnergy+2.0*electron_mass_c2);
1437 G4double cp1 = std::sqrt(cp1s);
1438 G4double cp = std::sqrt(cps);
1439 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1440
1441 //Distant longitudinal interactions
1442 G4double qm = 0.0;
1443
1444 if (resEnergy > kineticEnergy*(1e-6))
1445 {
1446 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1447 }
1448 else
1449 {
1450 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1451 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1452 }
1453
1454 if (qm < resEnergy)
1455 {
1456 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1457 }
1458 else
1459 {
1460 sdLong = 0.0;
1461 }
1462
1463 if (sdLong > 0) {
1464 sdTrans = std::max(std::log(gamma2)-beta2-densityCorrection,0.0);
1465 sdDist = sdTrans + sdLong;
1466 if (cutoffEnergy > resEnergy)
1467 {
1468 softCont = sdDist/resEnergy;
1469 }
1470 else
1471 {
1472 hardCont = sdDist/resEnergy;
1473 }
1474 }
1475
1476
1477 // Close collisions (Bhabha's cross section)
1478 G4double wl = std::max(cutoffEnergy,resEnergy);
1479 G4double wu = kineticEnergy;
1480
1481 if (wl < (wu-1*eV)) {
1482 hardCont += (1.0/wl)-(1.0/wu)-bha1*std::log(wu/wl)/kineticEnergy
1483 + bha2*(wu-wl)/(kineticEnergy*kineticEnergy)
1484 -bha3*((wu*wu)-(wl*wl))/(2.0*kineticEnergy*kineticEnergy*kineticEnergy)
1485 + bha4*((wu*wu*wu)-(wl*wl*wl))/(3.0*kineticEnergy*kineticEnergy*kineticEnergy*kineticEnergy);
1486 wu=wl;
1487 }
1488 wl = resEnergy;
1489 if (wl > (wu-1*eV))
1490 {
1491 theResult.first = softCont;
1492 theResult.second = hardCont;
1493 return theResult;
1494 }
1495 softCont += (1.0/wl)-(1.0/wu)-bha1*std::log(wu/wl)/kineticEnergy
1496 + bha2*(wu-wl)/(kineticEnergy*kineticEnergy)
1497 -bha3*((wu*wu)-(wl*wl))/(2.0*kineticEnergy*kineticEnergy*kineticEnergy)
1498 + bha4*((wu*wu*wu)-(wl*wl*wl))/(3.0*kineticEnergy*kineticEnergy*kineticEnergy*kineticEnergy);
1499
1500
1501 theResult.first = softCont;
1502 theResult.second = hardCont;
1503 return theResult;
1504
1505}
1506
1507G4double G4PenelopeIonisationModel::ComputeStoppingPowerForElectrons(G4double kinEnergy,
1508 G4double cutEnergy,
1509 G4double deltaFermi,
1510 G4double resEnergy)
1511{
1512 //Calculate constants
1513 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1514 G4double gamma2 = gamma*gamma;
1515 G4double beta2 = (gamma2-1.0)/gamma2;
1516 G4double cps = kinEnergy*(kinEnergy+2.0*electron_mass_c2);
1517 G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
1518 G4double sPower = 0.0;
1519 if (kinEnergy < resEnergy) return sPower;
1520
1521 //Distant interactions
1522 G4double cp1s = (kinEnergy-resEnergy)*(kinEnergy-resEnergy+2.0*electron_mass_c2);
1523 G4double cp1 = std::sqrt(cp1s);
1524 G4double cp = std::sqrt(cps);
1525 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1526
1527 //Distant longitudinal interactions
1528 G4double qm = 0.0;
1529
1530 if (resEnergy > kinEnergy*(1e-6))
1531 {
1532 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1533 }
1534 else
1535 {
1536 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1537 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1538 }
1539
1540 if (qm < resEnergy)
1541 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1542 else
1543 sdLong = 0.0;
1544
1545 if (sdLong > 0) {
1546 sdTrans = std::max(std::log(gamma2)-beta2-deltaFermi,0.0);
1547 sdDist = sdTrans + sdLong;
1548 if (cutEnergy > resEnergy) sPower = sdDist;
1549 }
1550
1551
1552 // Close collisions (Moeller's cross section)
1553 G4double wl = std::max(cutEnergy,resEnergy);
1554 G4double wu = 0.5*kinEnergy;
1555
1556 if (wl < (wu-1*eV)) wu=wl;
1557 wl = resEnergy;
1558 if (wl > (wu-1*eV)) return sPower;
1559 sPower += std::log(wu/wl)+(kinEnergy/(kinEnergy-wu))-(kinEnergy/(kinEnergy-wl))
1560 + (2.0 - amol)*std::log((kinEnergy-wu)/(kinEnergy-wl))
1561 + amol*((wu*wu)-(wl*wl))/(2.0*kinEnergy*kinEnergy);
1562 return sPower;
1563}
1564
1565
1566//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1567
1568G4double G4PenelopeIonisationModel::ComputeStoppingPowerForPositrons(G4double kinEnergy,
1569 G4double cutEnergy,
1570 G4double deltaFermi,
1571 G4double resEnergy)
1572{
1573 //Calculate constants
1574 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1575 G4double gamma2 = gamma*gamma;
1576 G4double beta2 = (gamma2-1.0)/gamma2;
1577 G4double cps = kinEnergy*(kinEnergy+2.0*electron_mass_c2);
1578 G4double amol = (kinEnergy/(kinEnergy+electron_mass_c2)) * (kinEnergy/(kinEnergy+electron_mass_c2));
1579 G4double help = (gamma+1.0)*(gamma+1.0);
1580 G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1581 G4double bha2 = amol*(3.0+1.0/help);
1582 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1583 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1584
1585 G4double sPower = 0.0;
1586 if (kinEnergy < resEnergy) return sPower;
1587
1588 //Distant interactions
1589 G4double cp1s = (kinEnergy-resEnergy)*(kinEnergy-resEnergy+2.0*electron_mass_c2);
1590 G4double cp1 = std::sqrt(cp1s);
1591 G4double cp = std::sqrt(cps);
1592 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1593
1594 //Distant longitudinal interactions
1595 G4double qm = 0.0;
1596
1597 if (resEnergy > kinEnergy*(1e-6))
1598 {
1599 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1600 }
1601 else
1602 {
1603 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1604 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1605 }
1606
1607 if (qm < resEnergy)
1608 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1609 else
1610 sdLong = 0.0;
1611
1612 if (sdLong > 0) {
1613 sdTrans = std::max(std::log(gamma2)-beta2-deltaFermi,0.0);
1614 sdDist = sdTrans + sdLong;
1615 if (cutEnergy > resEnergy) sPower = sdDist;
1616 }
1617
1618
1619 // Close collisions (Bhabha's cross section)
1620 G4double wl = std::max(cutEnergy,resEnergy);
1621 G4double wu = kinEnergy;
1622
1623 if (wl < (wu-1*eV)) wu=wl;
1624 wl = resEnergy;
1625 if (wl > (wu-1*eV)) return sPower;
1626 sPower += std::log(wu/wl)-bha1*(wu-wl)/kinEnergy
1627 + bha2*((wu*wu)-(wl*wl))/(2.0*kinEnergy*kinEnergy)
1628 - bha3*((wu*wu*wu)-(wl*wl*wl))/(3.0*kinEnergy*kinEnergy*kinEnergy)
1629 + bha4*((wu*wu*wu*wu)-(wl*wl*wl*wl))/(4.0*kinEnergy*kinEnergy*kinEnergy*kinEnergy);
1630 return sPower;
1631}
1632
1633//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1634
1635/* Notice: the methods here above are only temporary. They will become obsolete in a while */
1636
1637#include "G4VDataSetAlgorithm.hh"
1638#include "G4LinLogLogInterpolation.hh"
1639#include "G4CompositeEMDataSet.hh"
1640#include "G4EMDataSet.hh"
1641
1642std::vector<G4VEMDataSet*>* G4PenelopeIonisationModel::BuildCrossSectionTable(const
1643 G4ParticleDefinition* theParticle)
1644{
1645 std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;
1646
1647 size_t nOfBins = 200;
1648 G4PhysicsLogVector* theLogVector = new G4PhysicsLogVector(LowEnergyLimit(),
1649 HighEnergyLimit(),
1650 nOfBins);
1651 G4DataVector* energies;
1652 G4DataVector* cs;
1653
1654 const G4ProductionCutsTable* theCoupleTable=
1655 G4ProductionCutsTable::GetProductionCutsTable();
1656 size_t numOfCouples = theCoupleTable->GetTableSize();
1657
1658
1659 for (size_t m=0; m<numOfCouples; m++)
1660 {
1661 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
1662 const G4Material* material= couple->GetMaterial();
1663 const G4ElementVector* elementVector = material->GetElementVector();
1664 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
1665 G4double electronVolumeDensity =
1666 material->GetTotNbOfElectPerVolume(); //electron density
1667 G4int nElements = material->GetNumberOfElements();
1668
1669 G4double tcut = (*(theCoupleTable->GetEnergyCutsVector(1)))[couple->GetIndex()];
1670
1671 G4VDataSetAlgorithm* algo = new G4LinLogLogInterpolation();
1672
1673 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
1674
1675 for (G4int i=0; i<nElements; i++)
1676 {
1677 G4int iZ = (G4int) (*elementVector)[i]->GetZ();
1678 energies = new G4DataVector;
1679 cs = new G4DataVector;
1680 G4double density = nAtomsPerVolume[i];
1681 for (size_t bin=0; bin<nOfBins; bin++)
1682 {
1683 G4double e = theLogVector->GetLowEdgeEnergy(bin);
1684 energies->push_back(e);
1685 G4double value = 0.0;
1686 if(e > tcut)
1687 {
1688 G4double ratio = CalculateCrossSectionsRatio(e,tcut,iZ,
1689 electronVolumeDensity,
1690 theParticle);
1691 value = crossSectionHandler->FindValue(iZ,e)*ratio*density;
1692 }
1693 cs->push_back(value);
1694 }
1695 G4VDataSetAlgorithm* algo1 = algo->Clone();
1696 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,algo1,1.,1.);
1697 setForMat->AddComponent(elSet);
1698 }
1699 set->push_back(setForMat);
1700 }
1701 return set;
1702}
1703
1704//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1705
1706G4int G4PenelopeIonisationModel::SampleRandomAtom(const G4MaterialCutsCouple* couple,
1707 G4double e) const
1708{
1709 // Select randomly an element within the material, according to the weight
1710 // determined by the cross sections in the data set
1711
1712 const G4Material* material = couple->GetMaterial();
1713 G4int nElements = material->GetNumberOfElements();
1714
1715 // Special case: the material consists of one element
1716 if (nElements == 1)
1717 {
1718 G4int Z = (G4int) material->GetZ();
1719 return Z;
1720 }
1721
1722 // Composite material
1723 const G4ElementVector* elementVector = material->GetElementVector();
1724 size_t materialIndex = couple->GetIndex();
1725
1726 G4VEMDataSet* materialSet = (*theXSTable)[materialIndex];
1727 G4double materialCrossSection0 = 0.0;
1728 G4DataVector cross;
1729 cross.clear();
1730 for ( G4int i=0; i < nElements; i++ )
1731 {
1732 G4double cr = materialSet->GetComponent(i)->FindValue(e);
1733 materialCrossSection0 += cr;
1734 cross.push_back(materialCrossSection0);
1735 }
1736
1737 G4double random = G4UniformRand() * materialCrossSection0;
1738
1739 for (G4int k=0 ; k < nElements ; k++ )
1740 {
1741 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
1742 }
1743 // It should never get here
1744 return 0;
1745}
1746
1747
Note: See TracBrowser for help on using the repository browser.