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

Last change on this file since 1346 was 1340, checked in by garnier, 15 years ago

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