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

Last change on this file since 1058 was 1055, checked in by garnier, 17 years ago

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