source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermoreIonisationModel.cc

Last change on this file was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 25.6 KB
RevLine 
[1058]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//
[1347]26// $Id: G4LivermoreIonisationModel.cc,v 1.13 2010/12/02 16:06:29 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
[1058]28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 12 Jan 2009 L Pandola Migration from process to model
34// 03 Mar 2009 L Pandola Bug fix (release memory in the destructor)
35// 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
36// - apply internal high-energy limit only in constructor
37// - do not apply low-energy limit (default is 0)
38// - simplify sampling of deexcitation by using cut in energy
39// - set activation of Auger "false"
40// - remove initialisation of element selectors
41// 19 May 2009 L Pandola Explicitely set to zero pointers deleted in
42// Initialise(), since they might be checked later on
[1192]43// 23 Oct 2009 L Pandola
44// - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
45// set as "true" (default would be false)
[1340]46// 12 Oct 2010 L Pandola
47// - add debugging information about energy in
48// SampleDeexcitationAlongStep()
49// - generate fluorescence SampleDeexcitationAlongStep() only above
50// the cuts.
51//
[1058]52//
53
54#include "G4LivermoreIonisationModel.hh"
55#include "G4ParticleDefinition.hh"
56#include "G4MaterialCutsCouple.hh"
57#include "G4ProductionCutsTable.hh"
58#include "G4DynamicParticle.hh"
59#include "G4Element.hh"
60#include "G4AtomicTransitionManager.hh"
61#include "G4AtomicDeexcitation.hh"
62#include "G4AtomicShell.hh"
63#include "G4Gamma.hh"
64#include "G4Electron.hh"
65#include "G4CrossSectionHandler.hh"
66#include "G4ProcessManager.hh"
67#include "G4VEMDataSet.hh"
68#include "G4EMDataSet.hh"
69#include "G4eIonisationCrossSectionHandler.hh"
70#include "G4eIonisationSpectrum.hh"
71#include "G4VDataSetAlgorithm.hh"
72#include "G4SemiLogInterpolation.hh"
73#include "G4ShellVacancy.hh"
74#include "G4VDataSetAlgorithm.hh"
75#include "G4LogLogInterpolation.hh"
76#include "G4CompositeEMDataSet.hh"
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80
81G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
82 const G4String& nam)
83 :G4VEmModel(nam),isInitialised(false),crossSectionHandler(0),
84 energySpectrum(0),shellVacancy(0)
85{
86 fIntrinsicLowEnergyLimit = 10.0*eV;
87 fIntrinsicHighEnergyLimit = 100.0*GeV;
88 fNBinEnergyLoss = 360;
89 // SetLowEnergyLimit(fIntrinsicLowEnergyLimit);
90 SetHighEnergyLimit(fIntrinsicHighEnergyLimit);
91 //
92 verboseLevel = 0;
[1192]93 //By default: use deexcitation, not auger
94 SetDeexcitationFlag(true);
[1058]95 ActivateAuger(false);
96 //
[1315]97 //
[1058]98 // Notice: the fluorescence along step is generated only if it is
99 // set by the PROCESS (e.g. G4eIonisation) via the command
100 // process->ActivateDeexcitation(true);
101 //
102}
103
104//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
105
106G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
107{
108 if (energySpectrum) delete energySpectrum;
109 if (crossSectionHandler) delete crossSectionHandler;
110 if (shellVacancy) delete shellVacancy;
111}
112
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
114
115void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
116 const G4DataVector& cuts)
117{
118 //Check that the Livermore Ionisation is NOT attached to e+
119 if (particle != G4Electron::Electron())
120 {
121 G4cout << "ERROR: Livermore Ionisation Model is applicable only to electrons" << G4endl;
122 G4cout << "It cannot be registered to " << particle->GetParticleName() << G4endl;
123 G4Exception();
124 }
125
126 //Read energy spectrum
127 if (energySpectrum)
128 {
129 delete energySpectrum;
130 energySpectrum = 0;
131 }
132 energySpectrum = new G4eIonisationSpectrum();
[1340]133 if (verboseLevel > 3)
[1058]134 G4cout << "G4VEnergySpectrum is initialized" << G4endl;
135
136 //Initialize cross section handler
137 if (crossSectionHandler)
138 {
139 delete crossSectionHandler;
140 crossSectionHandler = 0;
141 }
142
143 G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
144 crossSectionHandler = new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
145 LowEnergyLimit(),HighEnergyLimit(),
146 fNBinEnergyLoss);
147 crossSectionHandler->Clear();
148 crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
149 //This is used to retrieve cross section values later on
[1347]150 G4VEMDataSet* emdata =
151 crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
152 //The method BuildMeanFreePathForMaterials() is required here only to force
153 //the building of an internal table: the output pointer can be deleted
154 delete emdata;
155
[1058]156 //Fluorescence data
157 transitionManager = G4AtomicTransitionManager::Instance();
158 if (shellVacancy) delete shellVacancy;
159 shellVacancy = new G4ShellVacancy();
160 InitialiseFluorescence();
161
162 if (verboseLevel > 0)
163 {
164 G4cout << "Livermore Ionisation model is initialized " << G4endl
165 << "Energy range: "
166 << LowEnergyLimit() / keV << " keV - "
167 << HighEnergyLimit() / GeV << " GeV"
168 << G4endl;
169 }
170
[1340]171 if (verboseLevel > 3)
[1058]172 {
173 G4cout << "Cross section data: " << G4endl;
174 crossSectionHandler->PrintData();
175 G4cout << "Parameters: " << G4endl;
176 energySpectrum->PrintData();
177 }
178
179 if(isInitialised) return;
180 fParticleChange = GetParticleChangeForLoss();
181 isInitialised = true;
182}
183
184//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
185
186G4double G4LivermoreIonisationModel::MinEnergyCut(const G4ParticleDefinition*,
187 const G4MaterialCutsCouple*)
188{
189 return 250.*eV;
190}
191
192//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
193
194G4double G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
195 G4double energy,
196 G4double Z, G4double,
197 G4double cutEnergy,
198 G4double)
199{
200 G4int iZ = (G4int) Z;
201 if (!crossSectionHandler)
202 {
203 G4cout << "G4LivermoreIonisationModel::ComputeCrossSectionPerAtom" << G4endl;
204 G4cout << "The cross section handler is not correctly initialized" << G4endl;
205 G4Exception();
[1347]206 return 0;
[1058]207 }
208
209 //The cut is already included in the crossSectionHandler
210 G4double cs =
211 crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,cutEnergy,iZ);
212
213 if (verboseLevel > 1)
214 {
215 G4cout << "G4LivermoreIonisationModel " << G4endl;
216 G4cout << "Cross section for delta emission > " << cutEnergy/keV << " keV at " <<
217 energy/keV << " keV and Z = " << iZ << " --> " << cs/barn << " barn" << G4endl;
218 }
219 return cs;
220}
221
222
223//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
224
225G4double G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
226 const G4ParticleDefinition* ,
227 G4double kineticEnergy,
228 G4double cutEnergy)
229{
230 G4double sPower = 0.0;
231
232 const G4ElementVector* theElementVector = material->GetElementVector();
233 size_t NumberOfElements = material->GetNumberOfElements() ;
234 const G4double* theAtomicNumDensityVector =
235 material->GetAtomicNumDensityVector();
236
237 // loop for elements in the material
238 for (size_t iel=0; iel<NumberOfElements; iel++ )
239 {
240 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
241 G4int nShells = transitionManager->NumberOfShells(iZ);
242 for (G4int n=0; n<nShells; n++)
243 {
244 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
245 kineticEnergy, n);
246 G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
247 sPower += e * cs * theAtomicNumDensityVector[iel];
248 }
249 G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
250 sPower += esp * theAtomicNumDensityVector[iel];
251 }
252
253 if (verboseLevel > 2)
254 {
255 G4cout << "G4LivermoreIonisationModel " << G4endl;
256 G4cout << "Stopping power < " << cutEnergy/keV << " keV at " <<
257 kineticEnergy/keV << " keV = " << sPower/(keV/mm) << " keV/mm" << G4endl;
258 }
259
260 return sPower;
261}
262
263//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
264
265void G4LivermoreIonisationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
266 const G4MaterialCutsCouple* couple,
267 const G4DynamicParticle* aDynamicParticle,
268 G4double cutE,
269 G4double maxE)
270{
271
272 G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
273
274 if (kineticEnergy <= fIntrinsicLowEnergyLimit)
275 {
276 fParticleChange->SetProposedKineticEnergy(0.);
277 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
278 return ;
279 }
280
281 // Select atom and shell
282 G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
283 G4int shell = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
284 const G4AtomicShell* atomicShell =
285 (G4AtomicTransitionManager::Instance())->Shell(Z, shell);
286 G4double bindingEnergy = atomicShell->BindingEnergy();
287 G4int shellId = atomicShell->ShellId();
288
289 // Sample delta energy using energy interval for delta-electrons
290 G4double energyMax =
291 std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
292 G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
293 kineticEnergy, shell);
294
295 if (energyDelta == 0.) //nothing happens
296 return;
297
298 // Transform to shell potential
299 G4double deltaKinE = energyDelta + 2.0*bindingEnergy;
300 G4double primaryKinE = kineticEnergy + 2.0*bindingEnergy;
301
302 // sampling of scattering angle neglecting atomic motion
303 G4double deltaMom = std::sqrt(deltaKinE*(deltaKinE + 2.0*electron_mass_c2));
304 G4double primaryMom = std::sqrt(primaryKinE*(primaryKinE + 2.0*electron_mass_c2));
305
306 G4double cost = deltaKinE * (primaryKinE + 2.0*electron_mass_c2)
307 / (deltaMom * primaryMom);
308 if (cost > 1.) cost = 1.;
309 G4double sint = std::sqrt(1. - cost*cost);
310 G4double phi = twopi * G4UniformRand();
311 G4double dirx = sint * std::cos(phi);
312 G4double diry = sint * std::sin(phi);
313 G4double dirz = cost;
314
315 // Rotate to incident electron direction
316 G4ThreeVector primaryDirection = aDynamicParticle->GetMomentumDirection();
317 G4ThreeVector deltaDir(dirx,diry,dirz);
318 deltaDir.rotateUz(primaryDirection);
319 //Updated components
320 dirx = deltaDir.x();
321 diry = deltaDir.y();
322 dirz = deltaDir.z();
323
324 // Take into account atomic motion del is relative momentum of the motion
325 // kinetic energy of the motion == bindingEnergy in V.Ivanchenko model
326 cost = 2.0*G4UniformRand() - 1.0;
327 sint = std::sqrt(1. - cost*cost);
328 phi = twopi * G4UniformRand();
329 G4double del = std::sqrt(bindingEnergy *(bindingEnergy + 2.0*electron_mass_c2))
330 / deltaMom;
331 dirx += del* sint * std::cos(phi);
332 diry += del* sint * std::sin(phi);
333 dirz += del* cost;
334
335 // Find out new primary electron direction
336 G4double finalPx = primaryMom*primaryDirection.x() - deltaMom*dirx;
337 G4double finalPy = primaryMom*primaryDirection.y() - deltaMom*diry;
338 G4double finalPz = primaryMom*primaryDirection.z() - deltaMom*dirz;
339
340 //Ok, ready to create the delta ray
341 G4DynamicParticle* theDeltaRay = new G4DynamicParticle();
342 theDeltaRay->SetKineticEnergy(energyDelta);
343 G4double norm = 1.0/std::sqrt(dirx*dirx + diry*diry + dirz*dirz);
344 dirx *= norm;
345 diry *= norm;
346 dirz *= norm;
347 theDeltaRay->SetMomentumDirection(dirx, diry, dirz);
348 theDeltaRay->SetDefinition(G4Electron::Electron());
349 fvect->push_back(theDeltaRay);
350
351 //This is the amount of energy available for fluorescence
352 G4double theEnergyDeposit = bindingEnergy;
353
354 // fill ParticleChange
355 // changed energy and momentum of the actual particle
356 G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
357 if(finalKinEnergy < 0.0)
358 {
359 theEnergyDeposit += finalKinEnergy;
360 finalKinEnergy = 0.0;
361 }
362 else
363 {
364 G4double norm = 1.0/std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz);
365 finalPx *= norm;
366 finalPy *= norm;
367 finalPz *= norm;
368 fParticleChange->ProposeMomentumDirection(finalPx, finalPy, finalPz);
369 }
370 fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
371
372 // deexcitation may be active per G4Region
373 if(DeexcitationFlag() && Z > 5)
374 {
375 G4ProductionCutsTable* theCoupleTable =
376 G4ProductionCutsTable::GetProductionCutsTable();
377 // Retrieve cuts for gammas
378 G4double cutG = (*theCoupleTable->GetEnergyCutsVector(0))[couple->GetIndex()];
379 if(theEnergyDeposit > cutG || theEnergyDeposit > cutE)
380 {
381 deexcitationManager.SetCutForSecondaryPhotons(cutG);
382 deexcitationManager.SetCutForAugerElectrons(cutE);
383 std::vector<G4DynamicParticle*>* secondaryVector =
384 deexcitationManager.GenerateParticles(Z, shellId);
385 G4DynamicParticle* aSecondary;
386
387 if (secondaryVector)
388 {
389 for (size_t i = 0; i<secondaryVector->size(); i++)
390 {
391 aSecondary = (*secondaryVector)[i];
392 //Check if it is a valid secondary
393 if (aSecondary)
394 {
395 G4double e = aSecondary->GetKineticEnergy();
396 if (e < theEnergyDeposit)
397 {
398 theEnergyDeposit -= e;
399 fvect->push_back(aSecondary);
[1347]400 aSecondary = 0;
401 (*secondaryVector)[i]=0;
[1058]402 }
403 else
404 {
405 delete aSecondary;
406 (*secondaryVector)[i] = 0;
407 }
408 }
409 }
[1347]410 //secondaryVector = 0;
411 delete secondaryVector;
[1058]412 }
413 }
414 }
415
416 if (theEnergyDeposit < 0)
417 {
418 G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
419 << theEnergyDeposit/eV << " eV" << G4endl;
420 theEnergyDeposit = 0.0;
421 }
422
423 //Assign local energy deposit
424 fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
425
426 if (verboseLevel > 1)
427 {
428 G4cout << "-----------------------------------------------------------" << G4endl;
429 G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
430 G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
431 G4cout << "-----------------------------------------------------------" << G4endl;
432 G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
433 G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
434 G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
435 G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
436 G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy+
437 theEnergyDeposit)/keV << " keV" << G4endl;
438 G4cout << "-----------------------------------------------------------" << G4endl;
439 }
440 return;
441}
442
443//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
444
445
446void G4LivermoreIonisationModel::SampleDeexcitationAlongStep(const G4Material* theMaterial,
447 const G4Track& theTrack,
448 G4double& eloss)
449{
450 //No call if there is no deexcitation along step
451 if (!DeexcitationFlag()) return;
452
453 //This method gets the energy loss calculated "Along the Step" and
454 //(including fluctuations) and produces explicit fluorescence/Auger
455 //secondaries. The eloss value is updated.
456 G4double energyLossBefore = eloss;
[1340]457
[1058]458 if (verboseLevel > 2)
[1340]459 {
460 G4cout << "-----------------------------------------------------------" << G4endl;
461 G4cout << " SampleDeexcitationAlongStep() from G4LivermoreIonisation" << G4endl;
462 G4cout << "Energy loss along step before deexcitation : " << energyLossBefore/keV <<
463 " keV" << G4endl;
464 }
[1058]465 G4double incidentEnergy = theTrack.GetDynamicParticle()->GetKineticEnergy();
466
467 G4ProductionCutsTable* theCoupleTable =
468 G4ProductionCutsTable::GetProductionCutsTable();
469 const G4MaterialCutsCouple* couple = theTrack.GetMaterialCutsCouple();
470 size_t index = couple->GetIndex();
471 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
472 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
473
474
475 std::vector<G4DynamicParticle*>* deexcitationProducts =
476 new std::vector<G4DynamicParticle*>;
477
478 if(eloss > cute || eloss > cutg)
479 {
480 const G4AtomicTransitionManager* transitionManager =
481 G4AtomicTransitionManager::Instance();
482 deexcitationManager.SetCutForSecondaryPhotons(cutg);
483 deexcitationManager.SetCutForAugerElectrons(cute);
484
485 size_t nElements = theMaterial->GetNumberOfElements();
486 const G4ElementVector* theElementVector = theMaterial->GetElementVector();
487
488 std::vector<G4DynamicParticle*>* secVector = 0;
489 G4DynamicParticle* aSecondary = 0;
490 //G4ParticleDefinition* type = 0;
491 G4double e;
492 G4ThreeVector position;
493 G4int shell, shellId;
494
495 // sample secondaries
496 G4double eTot = 0.0;
497 std::vector<G4int> n =
498 shellVacancy->GenerateNumberOfIonisations(couple,
499 incidentEnergy,eloss);
500 for (size_t i=0; i<nElements; i++)
501 {
502 G4int Z = (G4int)((*theElementVector)[i]->GetZ());
503 size_t nVacancies = n[i];
504 G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
505 if (nVacancies > 0 && Z > 5 && (maxE > cute || maxE > cutg))
506 {
507 for (size_t j=0; j<nVacancies; j++)
508 {
509 shell = crossSectionHandler->SelectRandomShell(Z, incidentEnergy);
510 shellId = transitionManager->Shell(Z, shell)->ShellId();
511 G4double maxEShell =
512 transitionManager->Shell(Z, shell)->BindingEnergy();
513 if (maxEShell > cute || maxEShell > cutg )
514 {
515 secVector = deexcitationManager.GenerateParticles(Z, shellId);
516 if (secVector)
517 {
518 for (size_t l = 0; l<secVector->size(); l++) {
519 aSecondary = (*secVector)[l];
520 if (aSecondary)
521 {
522 e = aSecondary->GetKineticEnergy();
[1340]523 G4double itsCut = cutg;
524 if (aSecondary->GetParticleDefinition() == G4Electron::Electron())
525 itsCut = cute;
526 if ( eTot + e <= eloss && e > itsCut )
[1058]527 {
528 eTot += e;
529 deexcitationProducts->push_back(aSecondary);
530 }
531 else
532 {
533 delete aSecondary;
534 }
535 }
536 }
537 delete secVector;
538 }
539 }
540 }
541 }
542 }
543 }
544
[1340]545 G4double energyLossInFluorescence = 0.0;
[1058]546 size_t nSecondaries = deexcitationProducts->size();
547 if (nSecondaries > 0)
548 {
[1340]549 //You may have already secondaries produced by SampleSubCutSecondaries()
550 //at the process G4VEnergyLossProcess
551 G4int secondariesBefore = fParticleChange->GetNumberOfSecondaries();
552 fParticleChange->SetNumberOfSecondaries(nSecondaries+secondariesBefore);
[1058]553 const G4StepPoint* preStep = theTrack.GetStep()->GetPreStepPoint();
554 const G4StepPoint* postStep = theTrack.GetStep()->GetPostStepPoint();
555 G4ThreeVector r = preStep->GetPosition();
556 G4ThreeVector deltaR = postStep->GetPosition();
557 deltaR -= r;
558 G4double t = preStep->GetGlobalTime();
559 G4double deltaT = postStep->GetGlobalTime();
560 deltaT -= t;
561 G4double time, q;
562 G4ThreeVector position;
[1340]563
[1058]564 for (size_t i=0; i<nSecondaries; i++)
565 {
566 G4DynamicParticle* part = (*deexcitationProducts)[i];
567 if (part)
568 {
569 G4double eSecondary = part->GetKineticEnergy();
570 eloss -= eSecondary;
571 if (eloss > 0.)
572 {
573 q = G4UniformRand();
574 time = deltaT*q + t;
575 position = deltaR*q;
576 position += r;
577 G4Track* newTrack = new G4Track(part, time, position);
[1340]578 energyLossInFluorescence += eSecondary;
[1058]579 pParticleChange->AddSecondary(newTrack);
580 }
581 else
582 {
583 eloss += eSecondary;
584 delete part;
585 part = 0;
586 }
587 }
588 }
589 }
590 delete deexcitationProducts;
591
[1340]592 //Check and verbosities. Ensure energy conservation
[1058]593 if (verboseLevel > 2)
[1340]594 {
595 G4cout << "Energy loss along step after deexcitation : " << eloss/keV <<
596 " keV" << G4endl;
597 }
598 if (verboseLevel > 1)
599 {
600 G4cout << "------------------------------------------------------------------" << G4endl;
601 G4cout << "Energy in fluorescence: " << energyLossInFluorescence/keV << " keV" << G4endl;
602 G4cout << "Residual energy loss: " << eloss/keV << " keV " << G4endl;
603 G4cout << "Total final: " << (energyLossInFluorescence+eloss)/keV << " keV" << G4endl;
604 G4cout << "Total initial: " << energyLossBefore/keV << " keV" << G4endl;
605 G4cout << "------------------------------------------------------------------" << G4endl;
606 }
607 if (verboseLevel > 0)
608 {
609 if (std::fabs(energyLossBefore-energyLossInFluorescence-eloss)>10*eV)
610 {
611 G4cout << "Found energy non-conservation at SampleDeexcitationAlongStep() " << G4endl;
612 G4cout << "Energy in fluorescence: " << energyLossInFluorescence/keV << " keV" << G4endl;
613 G4cout << "Residual energy loss: " << eloss/keV << " keV " << G4endl;
614 G4cout << "Total final: " << (energyLossInFluorescence+eloss)/keV << " keV" << G4endl;
615 G4cout << "Total initial: " << energyLossBefore/keV << " keV" << G4endl;
616 }
617 }
[1058]618}
619
620//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
621
622void G4LivermoreIonisationModel::InitialiseFluorescence()
623{
624 G4DataVector* ksi = 0;
625 G4DataVector* energyVector = 0;
626 size_t binForFluo = fNBinEnergyLoss/10;
627
[1347]628 //Used to produce a log-spaced energy grid. To be deleted at the end.
[1058]629 G4PhysicsLogVector* eVector = new G4PhysicsLogVector(LowEnergyLimit(),HighEnergyLimit(),
630 binForFluo);
631 const G4ProductionCutsTable* theCoupleTable=
632 G4ProductionCutsTable::GetProductionCutsTable();
633 size_t numOfCouples = theCoupleTable->GetTableSize();
634
635 // Loop on couples
636 for (size_t m=0; m<numOfCouples; m++)
637 {
638 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(m);
639 const G4Material* material= couple->GetMaterial();
640
641 const G4ElementVector* theElementVector = material->GetElementVector();
642 size_t NumberOfElements = material->GetNumberOfElements() ;
643 const G4double* theAtomicNumDensityVector =
644 material->GetAtomicNumDensityVector();
645
646 G4VDataSetAlgorithm* interp = new G4LogLogInterpolation();
647 G4VEMDataSet* xsis = new G4CompositeEMDataSet(interp, 1., 1.);
648 //loop on elements
649 G4double energyCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[m];
650 for (size_t iel=0; iel<NumberOfElements; iel++ )
651 {
652 G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
653 energyVector = new G4DataVector();
654 ksi = new G4DataVector();
655 //Loop on energy
656 for (size_t j = 0; j<binForFluo; j++)
657 {
658 G4double energy = eVector->GetLowEdgeEnergy(j);
659 G4double cross = 0.;
660 G4double eAverage= 0.;
661 G4int nShells = transitionManager->NumberOfShells(iZ);
662
663 for (G4int n=0; n<nShells; n++)
664 {
665 G4double e = energySpectrum->AverageEnergy(iZ, 0.0,energyCut,
666 energy, n);
667 G4double pro = energySpectrum->Probability(iZ, 0.0,energyCut,
668 energy, n);
669 G4double cs= crossSectionHandler->FindValue(iZ, energy, n);
670 eAverage += e * cs * theAtomicNumDensityVector[iel];
671 cross += cs * pro * theAtomicNumDensityVector[iel];
672 if(verboseLevel > 1)
673 {
674 G4cout << "Z= " << iZ
675 << " shell= " << n
676 << " E(keV)= " << energy/keV
677 << " Eav(keV)= " << e/keV
678 << " pro= " << pro
679 << " cs= " << cs
680 << G4endl;
681 }
682 }
683
684 G4double coeff = 0.0;
685 if(eAverage > 0.)
686 {
687 coeff = cross/eAverage;
688 eAverage /= cross;
689 }
690
691 if(verboseLevel > 1)
692 {
693 G4cout << "Ksi Coefficient for Z= " << iZ
694 << " E(keV)= " << energy/keV
695 << " Eav(keV)= " << eAverage/keV
696 << " coeff= " << coeff
697 << G4endl;
698 }
699 energyVector->push_back(energy);
700 ksi->push_back(coeff);
701 }
702 G4VEMDataSet* p = new G4EMDataSet(iZ,energyVector,ksi,interp,1.,1.);
703 xsis->AddComponent(p);
704 }
[1340]705 if(verboseLevel>3) xsis->PrintData();
[1058]706 shellVacancy->AddXsiTable(xsis);
707 }
[1347]708 if (eVector)
709 delete eVector;
[1058]710}
711
712//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
713
714void G4LivermoreIonisationModel::ActivateAuger(G4bool val)
715{
[1192]716 if (!DeexcitationFlag() && val)
717 {
718 G4cout << "WARNING - G4LivermoreIonisationModel" << G4endl;
719 G4cout << "The use of the Atomic Deexcitation Manager is set to false " << G4endl;
720 G4cout << "Therefore, Auger electrons will be not generated anyway" << G4endl;
721 }
[1058]722 deexcitationManager.ActivateAugerElectronProduction(val);
[1192]723 if (verboseLevel > 1)
724 G4cout << "Auger production set to " << val << G4endl;
[1058]725}
726
Note: See TracBrowser for help on using the repository browser.