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

Last change on this file since 973 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 60.6 KB
Line 
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: G4PenelopeIonisationModel.cc,v 1.3 2008/12/15 09:23:06 pandola Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
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;
519
520 G4double localEnergyDeposit = ionEnergy;
521 G4double energyInFluorescence = 0.0*eV;
522
523 std::vector<G4DynamicParticle*> *photonVector = 0;
524 if (fUseAtomicDeexcitation)
525 {
526 if (iZ>5 && (ionEnergy > cutG || ionEnergy > cutE))
527 {
528 photonVector = deexcitationManager.GenerateParticles(iZ,shellId);
529 //Check for single photons if they are above threshold
530 for (size_t k=0;k<photonVector->size();k++)
531 {
532 G4DynamicParticle* aPhoton = (*photonVector)[k];
533 if (aPhoton)
534 {
535 G4double itsCut = cutG;
536 if (aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cutE;
537 G4double itsEnergy = aPhoton->GetKineticEnergy();
538 if (itsEnergy > itsCut && itsEnergy <= ionEnergy)
539 {
540 localEnergyDeposit -= itsEnergy;
541 energyInFluorescence += itsEnergy;
542 }
543 else
544 {
545 delete aPhoton;
546 (*photonVector)[k] = 0;
547 }
548 }
549 }
550 }
551 }
552
553 // Generate the delta ray
554 G4double sin2 = std::sqrt(1. - cosThetaSecondary*cosThetaSecondary);
555 G4double phi2 = twopi * G4UniformRand();
556
557 G4double xEl = sin2 * std::cos(phi2);
558 G4double yEl = sin2 * std::sin(phi2);
559 G4double zEl = cosThetaSecondary;
560 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
561 eDirection.rotateUz(particleDirection0);
562
563 G4DynamicParticle* deltaElectron = new G4DynamicParticle (G4Electron::Electron(),
564 eDirection,eKineticEnergy) ;
565 fvect->push_back(deltaElectron);
566
567 //Generate fluorescence, if it is the case
568 //This block is executed only if there is at least one secondary photon produced by
569 //G4AtomicDeexcitation
570 if (photonVector)
571 {
572 for (size_t ll=0;ll<photonVector->size();ll++)
573 if ((*photonVector)[ll])
574 {
575 G4DynamicParticle* aFluorescencePhoton = (*photonVector)[ll];
576 fvect->push_back(aFluorescencePhoton);
577 }
578 }
579 delete photonVector;
580
581 if (localEnergyDeposit < 0)
582 {
583 G4cout << "WARNING-"
584 << "G4PenelopeIonisationModel::SampleSecondaries - Negative energy deposit"
585 << G4endl;
586 localEnergyDeposit=0.;
587 }
588 fParticleChange->ProposeLocalEnergyDeposit(localEnergyDeposit);
589
590 if (verboseLevel > 1)
591 {
592 G4cout << "-----------------------------------------------------------" << G4endl;
593 G4cout << "Energy balance from G4PenelopeIonisation" << G4endl;
594 G4cout << "Incoming primary energy: " << kineticEnergy0/keV << " keV" << G4endl;
595 G4cout << "-----------------------------------------------------------" << G4endl;
596 G4cout << "Outgoing primary energy: " << kineticEnergy1/keV << " keV" << G4endl;
597 G4cout << "Delta ray " << eKineticEnergy/keV << " keV" << G4endl;
598 G4cout << "Fluorescence: " << energyInFluorescence/keV << " keV" << G4endl;
599 G4cout << "Local energy deposit " << localEnergyDeposit/keV << " keV" << G4endl;
600 G4cout << "Total final state: " << (eKineticEnergy+energyInFluorescence+kineticEnergy1+
601 localEnergyDeposit)/keV <<
602 " keV" << G4endl;
603 G4cout << "-----------------------------------------------------------" << G4endl;
604 }
605 if (verboseLevel > 0)
606 {
607 G4double energyDiff = std::fabs(eKineticEnergy+energyInFluorescence+kineticEnergy1+
608 localEnergyDeposit-kineticEnergy0);
609 if (energyDiff > 0.05*keV)
610 G4cout << "Warning from G4PenelopeIonisation: problem with energy conservation: " <<
611 (eKineticEnergy+energyInFluorescence+kineticEnergy1+localEnergyDeposit)/keV <<
612 " keV (final) vs. " <<
613 kineticEnergy0/keV << " keV (initial)" << G4endl;
614 }
615}
616
617//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
618
619void G4PenelopeIonisationModel::ReadData()
620{
621 if (verboseLevel > 2)
622 {
623 G4cout << "Data from G4PenelopeIonisationModel read " << G4endl;
624 }
625 char* path = getenv("G4LEDATA");
626 if (!path)
627 {
628 G4String excep = "G4PenelopeIonisationModel - G4LEDATA environment variable not set!";
629 G4Exception(excep);
630 }
631 G4String pathString(path);
632 G4String pathFile = pathString + "/penelope/ion-pen.dat";
633 std::ifstream file(pathFile);
634
635 if (!file.is_open())
636 {
637 G4String excep = "G4PenelopeIonisationModel - data file " + pathFile + " not found!";
638 G4Exception(excep);
639 }
640
641 if (!ionizationEnergy || !resonanceEnergy || !occupationNumber || !shellFlag)
642 {
643 G4String excep = "G4PenelopeIonisationModel: problem with reading data from file";
644 G4Exception(excep);
645 }
646
647 G4int Z=1,nLevels=0;
648 G4int test,test1;
649
650 do{
651 G4DataVector* occVector = new G4DataVector;
652 G4DataVector* ionEVector = new G4DataVector;
653 G4DataVector* resEVector = new G4DataVector;
654 G4DataVector* shellIndVector = new G4DataVector;
655 file >> Z >> nLevels;
656 G4double a1,a2,a3,a4;
657 G4int k1,k2,k3;
658 for (G4int h=0;h<nLevels;h++)
659 {
660 //index,occup number,ion energy,res energy,fj0,kz,shell flag
661 file >> k1 >> a1 >> a2 >> a3 >> a4 >> k2 >> k3;
662 //Make explicit unit of measurements for ionisation and resonance energy,
663 // which is MeV
664 a2 *= MeV;
665 a3 *= MeV;
666 //
667 occVector->push_back(a1);
668 ionEVector->push_back(a2);
669 resEVector->push_back(a3);
670 shellIndVector->push_back((G4double) k3);
671 }
672 //Ok, done for element Z
673 occupationNumber->insert(std::make_pair(Z,occVector));
674 ionizationEnergy->insert(std::make_pair(Z,ionEVector));
675 resonanceEnergy->insert(std::make_pair(Z,resEVector));
676 shellFlag->insert(std::make_pair(Z,shellIndVector));
677 file >> test >> test1; //-1 -1 close the data for each Z
678 if (test > 0) {
679 G4String excep = "G4PenelopeIonisationModel - data file corrupted!";
680 G4Exception(excep);
681 }
682 }while (test != -2); //the very last Z is closed with -2 instead of -1
683}
684
685
686//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
687
688G4double G4PenelopeIonisationModel::CalculateDeltaFermi(G4double kinEnergy ,G4int Z,
689 G4double electronVolumeDensity)
690{
691 G4double plasmaEnergyCoefficient = 1.377e-39*(MeV*MeV*m3); //(e*hbar)^2/(epsilon0*electron_mass)
692 G4double plasmaEnergySquared = plasmaEnergyCoefficient*electronVolumeDensity;
693 // std::sqrt(plasmaEnergySquared) is the plasma energy of the solid (MeV)
694 G4double gam = 1.0+kinEnergy/electron_mass_c2;
695 G4double gam2=gam*gam;
696 G4double delta = 0.0;
697
698 //Density effect
699 G4double TST = ((G4double) Z)/(gam2*plasmaEnergySquared);
700
701 G4double wl2 = 0.0;
702 G4double fdel=0.0;
703 G4double wr=0;
704 size_t nbOsc = resonanceEnergy->find(Z)->second->size();
705 for(size_t i=0;i<nbOsc;i++)
706 {
707 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
708 wr = (*(resonanceEnergy->find(Z)->second))[i];
709 fdel += occupNb/(wr*wr+wl2);
710 }
711 if (fdel < TST) return delta;
712 G4double help1 = (*(resonanceEnergy->find(Z)->second))[nbOsc-1];
713 wl2 = help1*help1;
714 do{
715 wl2=wl2*2.0;
716 fdel = 0.0;
717 for (size_t ii=0;ii<nbOsc;ii++){
718 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[ii];
719 wr = (*(resonanceEnergy->find(Z)->second))[ii];
720 fdel += occupNb/(wr*wr+wl2);
721 }
722 }while (fdel > TST);
723 G4double wl2l=0.0;
724 G4double wl2u = wl2;
725 G4double control = 0.0;
726 do{
727 wl2=0.5*(wl2l+wl2u);
728 fdel = 0.0;
729 for (size_t jj=0;jj<nbOsc;jj++){
730 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[jj];
731 wr = (*(resonanceEnergy->find(Z)->second))[jj];
732 fdel += occupNb/(wr*wr+wl2);
733 }
734 if (fdel > TST)
735 wl2l = wl2;
736 else
737 wl2u = wl2;
738 control = wl2u-wl2l-wl2*1e-12;
739 }while(control>0);
740
741 //Density correction effect
742 for (size_t kk=0;kk<nbOsc;kk++){
743 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[kk];
744 wr = (*(resonanceEnergy->find(Z)->second))[kk];
745 delta += occupNb*std::log(1.0+wl2/(wr*wr));
746 }
747 delta = (delta/((G4double) Z))-wl2/(gam2*plasmaEnergySquared);
748 return delta;
749}
750
751//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
752
753void G4PenelopeIonisationModel::CalculateDiscreteForElectrons(G4double kinEnergy,G4double cutoffEnergy,
754 G4int Z,G4double electronVolumeDensity)
755{
756 if (verboseLevel > 2)
757 G4cout << "Entering in CalculateDiscreteForElectrons() for energy " <<
758 kinEnergy/keV << " keV " << G4endl;
759
760 //Initialization of variables to be calculated in this routine
761 kineticEnergy1=kinEnergy;
762 cosThetaPrimary=1.0;
763 energySecondary=0.0;
764 cosThetaSecondary=1.0;
765 iOsc=-1;
766
767 //constants
768 G4double rb=kinEnergy+2.0*electron_mass_c2;
769 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
770 G4double gamma2 = gamma*gamma;
771 G4double beta2 = (gamma2-1.0)/gamma2;
772 G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
773 G4double cps = kinEnergy*rb;
774 G4double cp = std::sqrt(cps);
775
776 G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
777 G4double distantTransvCS0 = std::max(std::log(gamma2)-beta2-delta,0.0);
778
779 G4double rl,rl1;
780
781 if (cutoffEnergy > kinEnergy) return; //delta rays are not generated
782
783 G4DataVector* qm = new G4DataVector();
784 G4DataVector* cumulHardCS = new G4DataVector();
785 std::vector<G4int> *typeOfInteraction = new std::vector<G4int>;
786 G4DataVector* nbOfLevel = new G4DataVector();
787
788 //Hard close collisions with outer shells
789 G4double wmaxc = 0.5*kinEnergy;
790 G4double closeCS0 = 0.0;
791 G4double closeCS = 0.0;
792 if (cutoffEnergy>0.1*eV)
793 {
794 rl=cutoffEnergy/kinEnergy;
795 rl1=1.0-rl;
796 if (rl < 0.5)
797 closeCS0 = (amol*(0.5-rl)+(1.0/rl)-(1.0/rl1)+(1.0-amol)*std::log(rl/rl1))/kinEnergy;
798 }
799
800 // Cross sections for the different oscillators
801
802 // totalHardCS contains the cumulative hard interaction cross section for the different
803 // excitable levels and the different interaction channels (close, distant, etc.),
804 // i.e.
805 // cumulHardCS[0] = 0.0
806 // cumulHardCS[1] = 1st excitable level (distant longitudinal only)
807 // cumulHardCS[2] = 1st excitable level (distant longitudinal + transverse)
808 // cumulHardCS[3] = 1st excitable level (distant longitudinal + transverse + close)
809 // cumulHardCS[4] = 1st excitable level (all channels) + 2nd excitable level (distant long only)
810 // etc.
811 // This is used for sampling the atomic level which is ionised and the channel of the
812 // interaction.
813 //
814 // For each index iFill of the cumulHardCS vector,
815 // nbOfLevel[iFill] contains the current excitable atomic level and
816 // typeOfInteraction[iFill] contains the current interaction channel, with the legenda:
817 // 1 = distant longitudinal interaction
818 // 2 = distant transverse interaction
819 // 3 = close collision
820 // 4 = close collision with outer shells (in this case nbOfLevel < 0 --> no binding energy)
821
822
823 G4int nOscil = ionizationEnergy->find(Z)->second->size();
824 G4double totalHardCS = 0.0;
825 G4double involvedElectrons = 0.0;
826 for (G4int i=0;i<nOscil;i++){
827 G4double wi = (*(resonanceEnergy->find(Z)->second))[i];
828 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
829
830 //Distant excitations
831 if (wi>cutoffEnergy && wi<kinEnergy)
832 {
833 if (wi>(1e-6*kinEnergy))
834 {
835 G4double cpp=std::sqrt((kinEnergy-wi)*(kinEnergy-wi+2.0*electron_mass_c2));
836 qm->push_back(std::sqrt((cp-cpp)*(cp-cpp)+electron_mass_c2*electron_mass_c2)-electron_mass_c2);
837 }
838 else
839 {
840 qm->push_back((wi*wi)/(beta2*2.0*electron_mass_c2));
841 }
842 if ((*qm)[i] < wi)
843 {
844
845 G4double distantLongitCS = occupNb*std::log(wi*((*qm)[i]+2.0*electron_mass_c2)/
846 ((*qm)[i]*(wi+2.0*electron_mass_c2)))/wi;
847 cumulHardCS->push_back(totalHardCS);
848 typeOfInteraction->push_back(1); //distant longitudinal
849 nbOfLevel->push_back((G4double) i); //only excitable level are counted
850 totalHardCS += distantLongitCS;
851
852 G4double distantTransvCS = occupNb*distantTransvCS0/wi;
853
854 cumulHardCS->push_back(totalHardCS);
855 typeOfInteraction->push_back(2); //distant tranverse
856 nbOfLevel->push_back((G4double) i);
857 totalHardCS += distantTransvCS;
858 }
859 }
860 else
861 qm->push_back(wi);
862
863 //close collisions
864 if(wi < wmaxc)
865 {
866 if (wi < cutoffEnergy)
867 involvedElectrons += occupNb;
868 else
869 {
870 rl=wi/kinEnergy;
871 rl1=1.0-rl;
872 closeCS = occupNb*(amol*(0.5-rl)+(1.0/rl)-(1.0/rl1)+(1.0-amol)*std::log(rl/rl1))/kinEnergy;
873 cumulHardCS->push_back(totalHardCS);
874 typeOfInteraction->push_back(3); //close
875 nbOfLevel->push_back((G4double) i);
876 totalHardCS += closeCS;
877 }
878 }
879 } // loop on the levels
880
881 cumulHardCS->push_back(totalHardCS);
882 typeOfInteraction->push_back(4); //close interaction with outer shells
883 nbOfLevel->push_back(-1.0);
884 totalHardCS += involvedElectrons*closeCS0;
885 cumulHardCS->push_back(totalHardCS); //this is the final value of the totalHardCS
886
887 if (totalHardCS < 1e-30) {
888 kineticEnergy1=kinEnergy;
889 cosThetaPrimary=1.0;
890 energySecondary=0.0;
891 cosThetaSecondary=0.0;
892 iOsc=-1;
893 delete qm;
894 delete cumulHardCS;
895 delete typeOfInteraction;
896 delete nbOfLevel;
897 return;
898 }
899
900 //Testing purposes
901 if (verboseLevel > 6)
902 {
903 for (size_t t=0;t<cumulHardCS->size();t++)
904 G4cout << (*cumulHardCS)[t] << " " << (*typeOfInteraction)[t] <<
905 " " << (*nbOfLevel)[t] << G4endl;
906 }
907
908 //Selection of the active oscillator on the basis of the cumulative cross sections
909 G4double TST = totalHardCS*G4UniformRand();
910 G4int is=0;
911 G4int js= nbOfLevel->size();
912 do{
913 G4int it=(is+js)/2;
914 if (TST > (*cumulHardCS)[it]) is=it;
915 if (TST <= (*cumulHardCS)[it]) js=it;
916 }while((js-is) > 1);
917
918 G4double UII=0.0;
919 G4double rkc=cutoffEnergy/kinEnergy;
920 G4double dde;
921 G4int kks;
922
923 G4int sampledInteraction = (*typeOfInteraction)[is];
924 iOsc = (G4int) (*nbOfLevel)[is];
925
926 if (verboseLevel > 4)
927 G4cout << "Chosen interaction #:" << sampledInteraction << " on level " << iOsc << G4endl;
928
929 //Generates the final state according to the sampled level and
930 //interaction channel
931
932 if (sampledInteraction == 1) //Hard distant longitudinal collisions
933 {
934 dde= (*(resonanceEnergy->find(Z)->second))[iOsc];
935 kineticEnergy1=kinEnergy-dde;
936 G4double qs=(*qm)[iOsc]/(1.0+((*qm)[iOsc]/(2.0*electron_mass_c2)));
937 G4double q=qs/(std::pow((qs/dde)*(1.0+(0.5*dde/electron_mass_c2)),G4UniformRand())-(0.5*qs/electron_mass_c2));
938 G4double qtrev = q*(q+2.0*electron_mass_c2);
939 G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
940 cosThetaPrimary = (cpps+cps-qtrev)/(2.0*cp*std::sqrt(cpps));
941 if (cosThetaPrimary>1.0) cosThetaPrimary=1.0;
942 //Energy and emission angle of the delta ray
943 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
944 //kks > 4 means that we are in an outer shell
945 if (kks>4)
946 energySecondary=dde;
947 else
948 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
949
950 cosThetaSecondary = 0.5*(dde*(kinEnergy+rb-dde)+qtrev)/std::sqrt(cps*qtrev);
951 if (cosThetaSecondary>1.0) cosThetaSecondary=1.0;
952 }
953 else if (sampledInteraction == 2) //Hard distant transverse collisions
954 {
955 dde=(*(resonanceEnergy->find(Z)->second))[iOsc];
956 kineticEnergy1=kinEnergy-dde;
957 cosThetaPrimary=1.0;
958 //Energy and emission angle of the delta ray
959 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
960 if (kks>4)
961 {
962 energySecondary=dde;
963 }
964 else
965 {
966 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
967 }
968 cosThetaSecondary = 1.0;
969 }
970 else if (sampledInteraction == 3 || sampledInteraction == 4) //Close interaction
971 {
972 if (sampledInteraction == 4) //interaction with inner shells
973 {
974 UII=0.0;
975 rkc = cutoffEnergy/kinEnergy;
976 iOsc = -1;
977 }
978 else
979 {
980 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
981 if (kks > 4) {
982 UII=0.0;
983 }
984 else
985 {
986 UII = (*(ionizationEnergy->find(Z)->second))[iOsc];
987 }
988 rkc = (*(resonanceEnergy->find(Z)->second))[iOsc]/kinEnergy;
989 }
990 G4double A = 0.5*amol;
991 G4double arkc = A*0.5*rkc;
992 G4double phi,rk2,rk,rkf;
993 do{
994 G4double fb = (1.0+arkc)*G4UniformRand();
995 if (fb<1.0)
996 {
997 rk=rkc/(1.0-fb*(1.0-(rkc*2.0)));
998 }
999 else{
1000 rk = rkc+(fb-1.0)*(0.5-rkc)/arkc;
1001 }
1002 rk2 = rk*rk;
1003 rkf = rk/(1.0-rk);
1004 phi = 1.0+(rkf*rkf)-rkf+amol*(rk2+rkf);
1005 }while ((G4UniformRand()*(1.0+A*rk2)) > phi);
1006 //Energy and scattering angle (primary electron);
1007 kineticEnergy1 = kinEnergy*(1.0-rk);
1008 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kinEnergy*(rb-(rk*kinEnergy))));
1009 //Energy and scattering angle of the delta ray
1010 energySecondary = kinEnergy-kineticEnergy1-UII;
1011 cosThetaSecondary = std::sqrt(rk*kinEnergy*rb/(kinEnergy*(rk*kinEnergy+2.0*electron_mass_c2)));
1012 }
1013 else
1014 {
1015 G4String excep = "G4PenelopeIonisationModel - Error in the calculation of the final state";
1016 G4Exception(excep);
1017 }
1018
1019 delete qm;
1020 delete cumulHardCS;
1021
1022 delete typeOfInteraction;
1023 delete nbOfLevel;
1024
1025 return;
1026}
1027
1028//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1029
1030 void G4PenelopeIonisationModel::CalculateDiscreteForPositrons(G4double kinEnergy,G4double cutoffEnergy,
1031 G4int Z,G4double electronVolumeDensity)
1032{
1033 kineticEnergy1=kinEnergy;
1034 cosThetaPrimary=1.0;
1035 energySecondary=0.0;
1036 cosThetaSecondary=1.0;
1037 iOsc=-1;
1038 //constants
1039 G4double rb=kinEnergy+2.0*electron_mass_c2;
1040 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1041 G4double gamma2 = gamma*gamma;
1042 G4double beta2 = (gamma2-1.0)/gamma2;
1043 G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
1044 G4double cps = kinEnergy*rb;
1045 G4double cp = std::sqrt(cps);
1046 G4double help = (gamma+1.0)*(gamma+1.0);
1047 G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1048 G4double bha2 = amol*(3.0+1.0/help);
1049 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1050 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1051
1052 G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
1053 G4double distantTransvCS0 = std::max(std::log(gamma2)-beta2-delta,0.0);
1054
1055 G4double rl,rl1;
1056
1057 if (cutoffEnergy > kinEnergy) return; //delta rays are not generated
1058
1059 G4DataVector* qm = new G4DataVector();
1060 G4DataVector* cumulHardCS = new G4DataVector();
1061 std::vector<G4int> *typeOfInteraction = new std::vector<G4int>;
1062 G4DataVector* nbOfLevel = new G4DataVector();
1063
1064
1065 //Hard close collisions with outer shells
1066 G4double wmaxc = kinEnergy;
1067 G4double closeCS0 = 0.0;
1068 G4double closeCS = 0.0;
1069 if (cutoffEnergy>0.1*eV)
1070 {
1071 rl=cutoffEnergy/kinEnergy;
1072 rl1=1.0-rl;
1073 if (rl < 1.0)
1074 closeCS0 = (((1.0/rl)-1.0) + bha1*std::log(rl) + bha2*rl1
1075 + (bha3/2.0)*((rl*rl)-1.0)
1076 + (bha4/3.0)*(1.0-(rl*rl*rl)))/kinEnergy;
1077 }
1078
1079 // Cross sections for the different oscillators
1080
1081 // totalHardCS contains the cumulative hard interaction cross section for the different
1082 // excitable levels and the different interaction channels (close, distant, etc.),
1083 // i.e.
1084 // cumulHardCS[0] = 0.0
1085 // cumulHardCS[1] = 1st excitable level (distant longitudinal only)
1086 // cumulHardCS[2] = 1st excitable level (distant longitudinal + transverse)
1087 // cumulHardCS[3] = 1st excitable level (distant longitudinal + transverse + close)
1088 // cumulHardCS[4] = 1st excitable level (all channels) + 2nd excitable level (distant long only)
1089 // etc.
1090 // This is used for sampling the atomic level which is ionised and the channel of the
1091 // interaction.
1092 //
1093 // For each index iFill of the cumulHardCS vector,
1094 // nbOfLevel[iFill] contains the current excitable atomic level and
1095 // typeOfInteraction[iFill] contains the current interaction channel, with the legenda:
1096 // 1 = distant longitudinal interaction
1097 // 2 = distant transverse interaction
1098 // 3 = close collision
1099 // 4 = close collision with outer shells (in this case nbOfLevel < 0 --> no binding energy)
1100
1101
1102 G4int nOscil = ionizationEnergy->find(Z)->second->size();
1103 G4double totalHardCS = 0.0;
1104 G4double involvedElectrons = 0.0;
1105 for (G4int i=0;i<nOscil;i++){
1106 G4double wi = (*(resonanceEnergy->find(Z)->second))[i];
1107 G4int occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
1108 //Distant excitations
1109 if (wi>cutoffEnergy && wi<kinEnergy)
1110 {
1111 if (wi>(1e-6*kinEnergy)){
1112 G4double cpp=std::sqrt((kinEnergy-wi)*(kinEnergy-wi+2.0*electron_mass_c2));
1113 qm->push_back(std::sqrt((cp-cpp)*(cp-cpp)+ electron_mass_c2 * electron_mass_c2)-electron_mass_c2);
1114 }
1115 else
1116 {
1117 qm->push_back(wi*wi/(beta2+2.0*electron_mass_c2));
1118 }
1119 //verificare che quando arriva qui il vettore ha SEMPRE l'i-esimo elemento
1120 if ((*qm)[i] < wi)
1121 {
1122
1123 G4double distantLongitCS = occupNb*std::log(wi*((*qm)[i]+2.0*electron_mass_c2)/
1124 ((*qm)[i]*(wi+2.0*electron_mass_c2)))/wi;
1125 cumulHardCS->push_back(totalHardCS);
1126 typeOfInteraction->push_back(1); //distant longitudinal
1127 nbOfLevel->push_back((G4double) i); //only excitable level are counted
1128 totalHardCS += distantLongitCS;
1129
1130 G4double distantTransvCS = occupNb*distantTransvCS0/wi;
1131
1132 cumulHardCS->push_back(totalHardCS);
1133 typeOfInteraction->push_back(2); //distant tranverse
1134 nbOfLevel->push_back((G4double) i);
1135 totalHardCS += distantTransvCS;
1136 }
1137 }
1138 else
1139 qm->push_back(wi);
1140
1141 //close collisions
1142 if(wi < wmaxc)
1143 {
1144 if (wi < cutoffEnergy) {
1145 involvedElectrons += occupNb;
1146 }
1147 else
1148 {
1149 rl=wi/kinEnergy;
1150 rl1=1.0-rl;
1151 closeCS = occupNb*(((1.0/rl)-1.0)+bha1*std::log(rl)+bha2*rl1
1152 + (bha3/2.0)*((rl*rl)-1.0)
1153 + (bha4/3.0)*(1.0-(rl*rl*rl)))/kinEnergy;
1154 cumulHardCS->push_back(totalHardCS);
1155 typeOfInteraction->push_back(3); //close
1156 nbOfLevel->push_back((G4double) i);
1157 totalHardCS += closeCS;
1158 }
1159 }
1160 } // loop on the levels
1161
1162 cumulHardCS->push_back(totalHardCS);
1163 typeOfInteraction->push_back(4); //close interaction with outer shells
1164 nbOfLevel->push_back(-1.0);
1165 totalHardCS += involvedElectrons*closeCS0;
1166 cumulHardCS->push_back(totalHardCS); //this is the final value of the totalHardCS
1167
1168 if (totalHardCS < 1e-30) {
1169 kineticEnergy1=kinEnergy;
1170 cosThetaPrimary=1.0;
1171 energySecondary=0.0;
1172 cosThetaSecondary=0.0;
1173 iOsc=-1;
1174 delete qm;
1175 delete cumulHardCS;
1176 delete typeOfInteraction;
1177 delete nbOfLevel;
1178 return;
1179 }
1180
1181 //Selection of the active oscillator on the basis of the cumulative cross sections
1182 G4double TST = totalHardCS*G4UniformRand();
1183 G4int is=0;
1184 G4int js= nbOfLevel->size();
1185 do{
1186 G4int it=(is+js)/2;
1187 if (TST > (*cumulHardCS)[it]) is=it;
1188 if (TST <= (*cumulHardCS)[it]) js=it;
1189 }while((js-is) > 1);
1190
1191 G4double UII=0.0;
1192 G4double rkc=cutoffEnergy/kinEnergy;
1193 G4double dde;
1194 G4int kks;
1195
1196 G4int sampledInteraction = (*typeOfInteraction)[is];
1197 iOsc = (G4int) (*nbOfLevel)[is];
1198
1199 //Generates the final state according to the sampled level and
1200 //interaction channel
1201
1202 if (sampledInteraction == 1) //Hard distant longitudinal collisions
1203 {
1204 dde= (*(resonanceEnergy->find(Z)->second))[iOsc];
1205 kineticEnergy1=kinEnergy-dde;
1206 G4double qs=(*qm)[iOsc]/(1.0+((*qm)[iOsc]/(2.0*electron_mass_c2)));
1207 G4double q=qs/(std::pow((qs/dde)*(1.0+(0.5*dde/electron_mass_c2)),G4UniformRand())-(0.5*qs/electron_mass_c2));
1208 G4double qtrev = q*(q+2.0*electron_mass_c2);
1209 G4double cpps = kineticEnergy1*(kineticEnergy1+2.0*electron_mass_c2);
1210 cosThetaPrimary = (cpps+cps-qtrev)/(2.0*cp*std::sqrt(cpps));
1211 if (cosThetaPrimary>1.0) cosThetaPrimary=1.0;
1212 //Energy and emission angle of the delta ray
1213 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1214 if (kks>4)
1215 energySecondary=dde;
1216
1217 else
1218 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
1219 cosThetaSecondary = 0.5*(dde*(kinEnergy+rb-dde)+qtrev)/std::sqrt(cps*qtrev);
1220 if (cosThetaSecondary>1.0) cosThetaSecondary=1.0;
1221 }
1222 else if (sampledInteraction == 2) //Hard distant transverse collisions
1223 {
1224 dde=(*(resonanceEnergy->find(Z)->second))[iOsc];
1225 kineticEnergy1=kinEnergy-dde;
1226 cosThetaPrimary=1.0;
1227 //Energy and emission angle of the delta ray
1228 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1229 if (kks>4)
1230 {
1231 energySecondary=dde;
1232 }
1233 else
1234 {
1235 energySecondary=dde-(*(ionizationEnergy->find(Z)->second))[iOsc];
1236 }
1237 cosThetaSecondary = 1.0;
1238 }
1239 else if (sampledInteraction == 3 || sampledInteraction == 4) //Close interaction
1240 {
1241 if (sampledInteraction == 4) //interaction with inner shells
1242 {
1243 UII=0.0;
1244 rkc = cutoffEnergy/kinEnergy;
1245 iOsc = -1;
1246 }
1247 else
1248 {
1249 kks = (G4int) (*(shellFlag->find(Z)->second))[iOsc];
1250 if (kks > 4) {
1251 UII=0.0;
1252 }
1253 else
1254 {
1255 UII = (*(ionizationEnergy->find(Z)->second))[iOsc];
1256 }
1257 rkc = (*(resonanceEnergy->find(Z)->second))[iOsc]/kinEnergy;
1258 }
1259 G4double phi,rk;
1260 do{
1261 rk=rkc/(1.0-G4UniformRand()*(1.0-rkc));
1262 phi = 1.0-rk*(bha1-rk*(bha2-rk*(bha3-bha4*rk)));
1263 }while ( G4UniformRand() > phi);
1264 //Energy and scattering angle (primary electron);
1265 kineticEnergy1 = kinEnergy*(1.0-rk);
1266 cosThetaPrimary = std::sqrt(kineticEnergy1*rb/(kinEnergy*(rb-(rk*kinEnergy))));
1267 //Energy and scattering angle of the delta ray
1268 energySecondary = kinEnergy-kineticEnergy1-UII;
1269 cosThetaSecondary = std::sqrt(rk*kinEnergy*rb/(kinEnergy*(rk*kinEnergy+2.0*electron_mass_c2)));
1270 }
1271 else
1272 {
1273 G4String excep = "G4PenelopeIonisationModel - Error in the calculation of the final state";
1274 G4Exception(excep);
1275 }
1276
1277 delete qm;
1278 delete cumulHardCS;
1279 delete typeOfInteraction;
1280 delete nbOfLevel;
1281
1282 return;
1283}
1284
1285//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1286
1287G4double G4PenelopeIonisationModel::CalculateCrossSectionsRatio(G4double kinEnergy,
1288 G4double cutoffEnergy,
1289 G4int Z, G4double electronVolumeDensity,
1290 const G4ParticleDefinition* theParticle)
1291{
1292 //Constants
1293 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1294 G4double gamma2 = gamma*gamma;
1295 G4double beta2 = (gamma2-1.0)/gamma2;
1296 G4double constant = pi*classic_electr_radius*classic_electr_radius*2.0*electron_mass_c2/beta2;
1297 G4double delta = CalculateDeltaFermi(kinEnergy,Z,electronVolumeDensity);
1298 G4int nbOsc = (G4int) resonanceEnergy->find(Z)->second->size();
1299 G4double softCS = 0.0;
1300 G4double hardCS = 0.0;
1301 for (G4int i=0;i<nbOsc;i++){
1302 G4double resEnergy = (*(resonanceEnergy->find(Z)->second))[i];
1303 std::pair<G4double,G4double> SoftAndHardXS(0.,0.);
1304 if (theParticle == G4Electron::Electron())
1305 SoftAndHardXS = CrossSectionsRatioForElectrons(kinEnergy,resEnergy,delta,cutoffEnergy);
1306 else if (theParticle == G4Positron::Positron())
1307 SoftAndHardXS = CrossSectionsRatioForPositrons(kinEnergy,resEnergy,delta,cutoffEnergy);
1308 G4double occupNb = (*(occupationNumber->find(Z)->second))[i];
1309 softCS += occupNb*constant*SoftAndHardXS.first;
1310 hardCS += occupNb*constant*SoftAndHardXS.second;
1311 }
1312 G4double ratio = 0.0;
1313
1314 if (softCS+hardCS) ratio = (hardCS)/(softCS+hardCS);
1315 return ratio;
1316}
1317
1318//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1319
1320std::pair<G4double,G4double> G4PenelopeIonisationModel::CrossSectionsRatioForElectrons(G4double kineticEnergy,
1321 G4double resEnergy,
1322 G4double densityCorrection,
1323 G4double cutoffEnergy)
1324{
1325 std::pair<G4double,G4double> theResult(0.,0.);
1326 if (kineticEnergy < resEnergy) return theResult;
1327
1328 //Calculate constants
1329 G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
1330 G4double gamma2 = gamma*gamma;
1331 G4double beta2 = (gamma2-1.0)/gamma2;
1332 G4double cps = kineticEnergy*(kineticEnergy+2.0*electron_mass_c2);
1333 G4double amol = (kineticEnergy/(kineticEnergy+electron_mass_c2)) * (kineticEnergy/(kineticEnergy+electron_mass_c2)) ;
1334 G4double hardCont = 0.0;
1335 G4double softCont = 0.0;
1336
1337 //Distant interactions
1338 G4double cp1s = (kineticEnergy-resEnergy)*(kineticEnergy-resEnergy+2.0*electron_mass_c2);
1339 G4double cp1 = std::sqrt(cp1s);
1340 G4double cp = std::sqrt(cps);
1341 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1342
1343 //Distant longitudinal interactions
1344 G4double qm = 0.0;
1345
1346 if (resEnergy > kineticEnergy*(1e-6))
1347 {
1348 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1349 }
1350 else
1351 {
1352 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1353 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1354 }
1355
1356 if (qm < resEnergy)
1357 {
1358 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1359 }
1360 else
1361 {
1362 sdLong = 0.0;
1363 }
1364
1365 if (sdLong > 0) {
1366 sdTrans = std::max(std::log(gamma2)-beta2-densityCorrection,0.0);
1367 sdDist = sdTrans + sdLong;
1368 if (cutoffEnergy > resEnergy)
1369 {
1370 softCont = sdDist/resEnergy;
1371 }
1372 else
1373 {
1374 hardCont = sdDist/resEnergy;
1375 }
1376 }
1377
1378 // Close collisions (Moeller's cross section)
1379 G4double wl = std::max(cutoffEnergy,resEnergy);
1380 G4double wu = 0.5*kineticEnergy;
1381
1382 if (wl < (wu-1*eV))
1383 {
1384 hardCont += (1.0/(kineticEnergy-wu))-(1.0/(kineticEnergy-wl))
1385 - (1.0/wu)+(1.0/wl)
1386 + (1.0-amol)*std::log(((kineticEnergy-wu)*wl)/((kineticEnergy-wl)*wu))/kineticEnergy
1387 + amol*(wu-wl)/(kineticEnergy*kineticEnergy);
1388 wu=wl;
1389 }
1390
1391 wl = resEnergy;
1392 if (wl > (wu-1*eV))
1393 {
1394 theResult.first = softCont;
1395 theResult.second = hardCont;
1396 return theResult;
1397 }
1398 softCont += (1.0/(kineticEnergy-wu))-(1.0/(kineticEnergy-wl))
1399 - (1.0/wu)+(1.0/wl)
1400 + (1.0-amol)*std::log(((kineticEnergy-wu)*wl)/((kineticEnergy-wl)*wu))/kineticEnergy
1401 + amol*(wu-wl)/(kineticEnergy*kineticEnergy);
1402 theResult.first = softCont;
1403 theResult.second = hardCont;
1404 return theResult;
1405}
1406
1407
1408//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1409
1410std::pair<G4double,G4double> G4PenelopeIonisationModel::CrossSectionsRatioForPositrons(G4double kineticEnergy,
1411 G4double resEnergy,
1412 G4double densityCorrection,
1413 G4double cutoffEnergy)
1414{
1415
1416 std::pair<G4double,G4double> theResult(0.,0.);
1417 if (kineticEnergy < resEnergy) return theResult;
1418
1419 //Calculate constants
1420 G4double gamma = 1.0+kineticEnergy/electron_mass_c2;
1421 G4double gamma2 = gamma*gamma;
1422 G4double beta2 = (gamma2-1.0)/gamma2;
1423 G4double cps = kineticEnergy*(kineticEnergy+2.0*electron_mass_c2);
1424 G4double amol = (kineticEnergy/(kineticEnergy+electron_mass_c2)) * (kineticEnergy/(kineticEnergy+electron_mass_c2)) ;
1425 G4double help = (gamma+1.0)*(gamma+1.0);
1426 G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1427 G4double bha2 = amol*(3.0+1.0/help);
1428 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1429 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1430 G4double hardCont = 0.0;
1431 G4double softCont = 0.0;
1432
1433 //Distant interactions
1434 G4double cp1s = (kineticEnergy-resEnergy)*(kineticEnergy-resEnergy+2.0*electron_mass_c2);
1435 G4double cp1 = std::sqrt(cp1s);
1436 G4double cp = std::sqrt(cps);
1437 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1438
1439 //Distant longitudinal interactions
1440 G4double qm = 0.0;
1441
1442 if (resEnergy > kineticEnergy*(1e-6))
1443 {
1444 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1445 }
1446 else
1447 {
1448 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1449 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1450 }
1451
1452 if (qm < resEnergy)
1453 {
1454 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1455 }
1456 else
1457 {
1458 sdLong = 0.0;
1459 }
1460
1461 if (sdLong > 0) {
1462 sdTrans = std::max(std::log(gamma2)-beta2-densityCorrection,0.0);
1463 sdDist = sdTrans + sdLong;
1464 if (cutoffEnergy > resEnergy)
1465 {
1466 softCont = sdDist/resEnergy;
1467 }
1468 else
1469 {
1470 hardCont = sdDist/resEnergy;
1471 }
1472 }
1473
1474
1475 // Close collisions (Bhabha's cross section)
1476 G4double wl = std::max(cutoffEnergy,resEnergy);
1477 G4double wu = kineticEnergy;
1478
1479 if (wl < (wu-1*eV)) {
1480 hardCont += (1.0/wl)-(1.0/wu)-bha1*std::log(wu/wl)/kineticEnergy
1481 + bha2*(wu-wl)/(kineticEnergy*kineticEnergy)
1482 -bha3*((wu*wu)-(wl*wl))/(2.0*kineticEnergy*kineticEnergy*kineticEnergy)
1483 + bha4*((wu*wu*wu)-(wl*wl*wl))/(3.0*kineticEnergy*kineticEnergy*kineticEnergy*kineticEnergy);
1484 wu=wl;
1485 }
1486 wl = resEnergy;
1487 if (wl > (wu-1*eV))
1488 {
1489 theResult.first = softCont;
1490 theResult.second = hardCont;
1491 return theResult;
1492 }
1493 softCont += (1.0/wl)-(1.0/wu)-bha1*std::log(wu/wl)/kineticEnergy
1494 + bha2*(wu-wl)/(kineticEnergy*kineticEnergy)
1495 -bha3*((wu*wu)-(wl*wl))/(2.0*kineticEnergy*kineticEnergy*kineticEnergy)
1496 + bha4*((wu*wu*wu)-(wl*wl*wl))/(3.0*kineticEnergy*kineticEnergy*kineticEnergy*kineticEnergy);
1497
1498
1499 theResult.first = softCont;
1500 theResult.second = hardCont;
1501 return theResult;
1502
1503}
1504
1505G4double G4PenelopeIonisationModel::ComputeStoppingPowerForElectrons(G4double kinEnergy,
1506 G4double cutEnergy,
1507 G4double deltaFermi,
1508 G4double resEnergy)
1509{
1510 //Calculate constants
1511 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1512 G4double gamma2 = gamma*gamma;
1513 G4double beta2 = (gamma2-1.0)/gamma2;
1514 G4double cps = kinEnergy*(kinEnergy+2.0*electron_mass_c2);
1515 G4double amol = (gamma-1.0)*(gamma-1.0)/gamma2;
1516 G4double sPower = 0.0;
1517 if (kinEnergy < resEnergy) return sPower;
1518
1519 //Distant interactions
1520 G4double cp1s = (kinEnergy-resEnergy)*(kinEnergy-resEnergy+2.0*electron_mass_c2);
1521 G4double cp1 = std::sqrt(cp1s);
1522 G4double cp = std::sqrt(cps);
1523 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1524
1525 //Distant longitudinal interactions
1526 G4double qm = 0.0;
1527
1528 if (resEnergy > kinEnergy*(1e-6))
1529 {
1530 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1531 }
1532 else
1533 {
1534 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1535 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1536 }
1537
1538 if (qm < resEnergy)
1539 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1540 else
1541 sdLong = 0.0;
1542
1543 if (sdLong > 0) {
1544 sdTrans = std::max(std::log(gamma2)-beta2-deltaFermi,0.0);
1545 sdDist = sdTrans + sdLong;
1546 if (cutEnergy > resEnergy) sPower = sdDist;
1547 }
1548
1549
1550 // Close collisions (Moeller's cross section)
1551 G4double wl = std::max(cutEnergy,resEnergy);
1552 G4double wu = 0.5*kinEnergy;
1553
1554 if (wl < (wu-1*eV)) wu=wl;
1555 wl = resEnergy;
1556 if (wl > (wu-1*eV)) return sPower;
1557 sPower += std::log(wu/wl)+(kinEnergy/(kinEnergy-wu))-(kinEnergy/(kinEnergy-wl))
1558 + (2.0 - amol)*std::log((kinEnergy-wu)/(kinEnergy-wl))
1559 + amol*((wu*wu)-(wl*wl))/(2.0*kinEnergy*kinEnergy);
1560 return sPower;
1561}
1562
1563
1564//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1565
1566G4double G4PenelopeIonisationModel::ComputeStoppingPowerForPositrons(G4double kinEnergy,
1567 G4double cutEnergy,
1568 G4double deltaFermi,
1569 G4double resEnergy)
1570{
1571 //Calculate constants
1572 G4double gamma = 1.0+kinEnergy/electron_mass_c2;
1573 G4double gamma2 = gamma*gamma;
1574 G4double beta2 = (gamma2-1.0)/gamma2;
1575 G4double cps = kinEnergy*(kinEnergy+2.0*electron_mass_c2);
1576 G4double amol = (kinEnergy/(kinEnergy+electron_mass_c2)) * (kinEnergy/(kinEnergy+electron_mass_c2));
1577 G4double help = (gamma+1.0)*(gamma+1.0);
1578 G4double bha1 = amol*(2.0*help-1.0)/(gamma2-1.0);
1579 G4double bha2 = amol*(3.0+1.0/help);
1580 G4double bha3 = amol*2.0*gamma*(gamma-1.0)/help;
1581 G4double bha4 = amol*(gamma-1.0)*(gamma-1.0)/help;
1582
1583 G4double sPower = 0.0;
1584 if (kinEnergy < resEnergy) return sPower;
1585
1586 //Distant interactions
1587 G4double cp1s = (kinEnergy-resEnergy)*(kinEnergy-resEnergy+2.0*electron_mass_c2);
1588 G4double cp1 = std::sqrt(cp1s);
1589 G4double cp = std::sqrt(cps);
1590 G4double sdLong=0.0, sdTrans = 0.0, sdDist=0.0;
1591
1592 //Distant longitudinal interactions
1593 G4double qm = 0.0;
1594
1595 if (resEnergy > kinEnergy*(1e-6))
1596 {
1597 qm = std::sqrt((cp-cp1)*(cp-cp1)+(electron_mass_c2*electron_mass_c2))-electron_mass_c2;
1598 }
1599 else
1600 {
1601 qm = resEnergy*resEnergy/(beta2*2.0*electron_mass_c2);
1602 qm = qm*(1.0-0.5*qm/electron_mass_c2);
1603 }
1604
1605 if (qm < resEnergy)
1606 sdLong = std::log(resEnergy*(qm+2.0*electron_mass_c2)/(qm*(resEnergy+2.0*electron_mass_c2)));
1607 else
1608 sdLong = 0.0;
1609
1610 if (sdLong > 0) {
1611 sdTrans = std::max(std::log(gamma2)-beta2-deltaFermi,0.0);
1612 sdDist = sdTrans + sdLong;
1613 if (cutEnergy > resEnergy) sPower = sdDist;
1614 }
1615
1616
1617 // Close collisions (Bhabha's cross section)
1618 G4double wl = std::max(cutEnergy,resEnergy);
1619 G4double wu = kinEnergy;
1620
1621 if (wl < (wu-1*eV)) wu=wl;
1622 wl = resEnergy;
1623 if (wl > (wu-1*eV)) return sPower;
1624 sPower += std::log(wu/wl)-bha1*(wu-wl)/kinEnergy
1625 + bha2*((wu*wu)-(wl*wl))/(2.0*kinEnergy*kinEnergy)
1626 - bha3*((wu*wu*wu)-(wl*wl*wl))/(3.0*kinEnergy*kinEnergy*kinEnergy)
1627 + bha4*((wu*wu*wu*wu)-(wl*wl*wl*wl))/(4.0*kinEnergy*kinEnergy*kinEnergy*kinEnergy);
1628 return sPower;
1629}
1630
1631//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1632
1633/* Notice: the methods here above are only temporary. They will become obsolete in a while */
1634
1635#include "G4VDataSetAlgorithm.hh"
1636#include "G4LinLogLogInterpolation.hh"
1637#include "G4CompositeEMDataSet.hh"
1638#include "G4EMDataSet.hh"
1639
1640std::vector<G4VEMDataSet*>* G4PenelopeIonisationModel::BuildCrossSectionTable(const
1641 G4ParticleDefinition* theParticle)
1642{
1643 std::vector<G4VEMDataSet*>* set = new std::vector<G4VEMDataSet*>;
1644
1645 size_t nOfBins = 200;
1646 G4PhysicsLogVector* theLogVector = new G4PhysicsLogVector(LowEnergyLimit(),
1647 HighEnergyLimit(),
1648 nOfBins);
1649 G4DataVector* energies;
1650 G4DataVector* cs;
1651
1652 const G4ProductionCutsTable* theCoupleTable=
1653 G4ProductionCutsTable::GetProductionCutsTable();
1654 size_t numOfCouples = theCoupleTable->GetTableSize();
1655
1656
1657 for (size_t m=0; m<numOfCouples; m++)
1658 {
1659 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
1660 const G4Material* material= couple->GetMaterial();
1661 const G4ElementVector* elementVector = material->GetElementVector();
1662 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
1663 G4double electronVolumeDensity =
1664 material->GetTotNbOfElectPerVolume(); //electron density
1665 G4int nElements = material->GetNumberOfElements();
1666
1667 G4double tcut = (*(theCoupleTable->GetEnergyCutsVector(1)))[couple->GetIndex()];
1668
1669 G4VDataSetAlgorithm* algo = new G4LinLogLogInterpolation();
1670
1671 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
1672
1673 for (G4int i=0; i<nElements; i++)
1674 {
1675 G4int iZ = (G4int) (*elementVector)[i]->GetZ();
1676 energies = new G4DataVector;
1677 cs = new G4DataVector;
1678 G4double density = nAtomsPerVolume[i];
1679 for (size_t bin=0; bin<nOfBins; bin++)
1680 {
1681 G4double e = theLogVector->GetLowEdgeEnergy(bin);
1682 energies->push_back(e);
1683 G4double value = 0.0;
1684 if(e > tcut)
1685 {
1686 G4double ratio = CalculateCrossSectionsRatio(e,tcut,iZ,
1687 electronVolumeDensity,
1688 theParticle);
1689 value = crossSectionHandler->FindValue(iZ,e)*ratio*density;
1690 }
1691 cs->push_back(value);
1692 }
1693 G4VDataSetAlgorithm* algo1 = algo->Clone();
1694 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,cs,algo1,1.,1.);
1695 setForMat->AddComponent(elSet);
1696 }
1697 set->push_back(setForMat);
1698 }
1699 return set;
1700}
1701
1702//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...
1703
1704G4int G4PenelopeIonisationModel::SampleRandomAtom(const G4MaterialCutsCouple* couple,
1705 G4double e) const
1706{
1707 // Select randomly an element within the material, according to the weight
1708 // determined by the cross sections in the data set
1709
1710 const G4Material* material = couple->GetMaterial();
1711 G4int nElements = material->GetNumberOfElements();
1712
1713 // Special case: the material consists of one element
1714 if (nElements == 1)
1715 {
1716 G4int Z = (G4int) material->GetZ();
1717 return Z;
1718 }
1719
1720 // Composite material
1721 const G4ElementVector* elementVector = material->GetElementVector();
1722 size_t materialIndex = couple->GetIndex();
1723
1724 G4VEMDataSet* materialSet = (*theXSTable)[materialIndex];
1725 G4double materialCrossSection0 = 0.0;
1726 G4DataVector cross;
1727 cross.clear();
1728 for ( G4int i=0; i < nElements; i++ )
1729 {
1730 G4double cr = materialSet->GetComponent(i)->FindValue(e);
1731 materialCrossSection0 += cr;
1732 cross.push_back(materialCrossSection0);
1733 }
1734
1735 G4double random = G4UniformRand() * materialCrossSection0;
1736
1737 for (G4int k=0 ; k < nElements ; k++ )
1738 {
1739 if (random <= cross[k]) return (G4int) (*elementVector)[k]->GetZ();
1740 }
1741 // It should never get here
1742 return 0;
1743}
1744
1745
Note: See TracBrowser for help on using the repository browser.