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

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

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