| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // File name: G4LowEnergyIonisation
30 //
31 // Author: Alessandra Forti, Vladimir I
32 //
33 // Creation date: March 1999
34 //
35 // Modifications:
36 // - 11.04.2000 VL
37 // Changing use of float and G4float casts t
38 // because of problems with optimisation (bu
39 // 10.04.2000 VL
40 // - Correcting Fluorescence transition probab
41 // non-radiative transitions. No Auger elect
42 // 10.04.2000 VL
43 // - Correction of incident electron final mom
44 // 07.04.2000 VL+LU
45 // - First implementation of continuous energy
46 // 22.03.2000 VL
47 // - 1 bug corrected in SelectRandomAtom metho
48 // 17.02.2000 Veronique Lefebure
49 // - 5 bugs corrected:
50 // *in Fluorescence, 2 bugs affecting
51 // . localEnergyDeposition and
52 // . number of emitted photons that was then
53 // *in EnergySampling method:
54 // . expon = Parms[13]+1; (instead of uncorr
55 // . rejection /= Parms[6];(instead of uncor
56 // . Parms[6] is apparently corrupted in the
57 // -->Compute normalisation into local var
58 // and use rejectionMax in stead of Parms
59 //
60 // Added Livermore data table construction met
61 // Modified BuildMeanFreePath to read new data
62 // Added EnergySampling method A. Forti
63 // Modified PostStepDoIt to insert sampling wi
64 // Added SelectRandomAtom A. Forti
65 // Added map of the elements A. Forti
66 // 20.09.00 V.Ivanchenko update fluctuations
67 // 24.04.01 V.Ivanchenko remove RogueWave
68 // 22.05.01 V.Ivanchenko update calculation of
69 // clean up the code
70 // 02.08.01 V.Ivanchenko fix energy conservati
71 // 18.08.01 V.Ivanchenko fix energy conservati
72 // 01.10.01 E. Guardincerri Replaced fluoresce
73 // according to desig
74 // 04.10.01 MGP Minor clean-up in
75 // compilation warnin
76 // prevent from acces
77 // 29.09.01 V.Ivanchenko revision based on
78 // 10.10.01 MGP Revision to improv
79 // consistency with d
80 // 18.10.01 V.Ivanchenko Add fluorescence A
81 // 18.10.01 MGP Revision to improv
82 // consistency with d
83 // 19.10.01 V.Ivanchenko update according t
84 // 26.10.01 V.Ivanchenko clean up deexcitat
85 // 28.10.01 V.Ivanchenko update printout
86 // 29.11.01 V.Ivanchenko New parametrisatio
87 // 25.03.02 V.Ivanchneko Fix in fluorescenc
88 // 28.03.02 V.Ivanchenko Add flag of fluore
89 // 28.05.02 V.Ivanchenko Remove flag fStopA
90 // 31.05.02 V.Ivanchenko Add path of Fluo +
91 // AtomicDeexcitation
92 // 03.06.02 MGP Restore fStopAndKi
93 // 19.06.02 VI Additional printou
94 // 30.07.02 VI Fix in restricted
95 // 20.09.02 VI Remove ActivateFlu
96 // 21.01.03 VI Cut per region
97 // 12.02.03 VI Change signature f
98 // 12.04.03 V.Ivanchenko Cut per region for
99 // 31.08.04 V.Ivanchenko Add density correc
100 //
101 // -------------------------------------------
102
103 #include "G4LowEnergyIonisation.hh"
104 #include "G4PhysicalConstants.hh"
105 #include "G4SystemOfUnits.hh"
106 #include "G4RDeIonisationSpectrum.hh"
107 #include "G4RDeIonisationCrossSectionHandler.h
108 #include "G4RDAtomicTransitionManager.hh"
109 #include "G4RDAtomicShell.hh"
110 #include "G4RDVDataSetAlgorithm.hh"
111 #include "G4RDSemiLogInterpolation.hh"
112 #include "G4RDLogLogInterpolation.hh"
113 #include "G4RDEMDataSet.hh"
114 #include "G4RDVEMDataSet.hh"
115 #include "G4RDCompositeEMDataSet.hh"
116 #include "G4EnergyLossTables.hh"
117 #include "G4RDShellVacancy.hh"
118 #include "G4UnitsTable.hh"
119 #include "G4Electron.hh"
120 #include "G4Gamma.hh"
121 #include "G4ProductionCutsTable.hh"
122
123 G4LowEnergyIonisation::G4LowEnergyIonisation(c
124 : G4eLowEnergyLoss(nam),
125 crossSectionHandler(0),
126 theMeanFreePath(0),
127 energySpectrum(0),
128 shellVacancy(0)
129 {
130 cutForPhotons = 250.0*eV;
131 cutForElectrons = 250.0*eV;
132 verboseLevel = 0;
133 }
134
135
136 G4LowEnergyIonisation::~G4LowEnergyIonisation(
137 {
138 delete crossSectionHandler;
139 delete energySpectrum;
140 delete theMeanFreePath;
141 delete shellVacancy;
142 }
143
144
145 void G4LowEnergyIonisation::BuildPhysicsTable(
146 {
147 if(verboseLevel > 0) {
148 G4cout << "G4LowEnergyIonisation::BuildPhy
149 << G4endl;
150 }
151
152 cutForDelta.clear();
153
154 // Create and fill IonisationParameters once
155 if( energySpectrum != 0 ) delete energySpect
156 energySpectrum = new G4RDeIonisationSpectrum
157
158 if(verboseLevel > 0) {
159 G4cout << "G4RDVEnergySpectrum is initiali
160 << G4endl;
161 }
162
163 // Create and fill G4RDCrossSectionHandler o
164
165 if ( crossSectionHandler != 0 ) delete cross
166 G4RDVDataSetAlgorithm* interpolation = new G
167 G4double lowKineticEnergy = GetLowerBoundEl
168 G4double highKineticEnergy = GetUpperBoundEl
169 G4int totBin = GetNbinEloss();
170 crossSectionHandler = new G4RDeIonisationCro
171 interpolation,
172 lowKineticEnergy,
173 highKineticEnergy,
174 totBin);
175 crossSectionHandler->LoadShellData("ioni/ion
176
177 if (verboseLevel > 0) {
178 G4cout << GetProcessName()
179 << " is created; Cross section data
180 << G4endl;
181 crossSectionHandler->PrintData();
182 G4cout << "Parameters: "
183 << G4endl;
184 energySpectrum->PrintData();
185 }
186
187 // Build loss table for IonisationIV
188
189 BuildLossTable(aParticleType);
190
191 if(verboseLevel > 0) {
192 G4cout << "The loss table is built"
193 << G4endl;
194 }
195
196 if (&aParticleType==G4Electron::Electron())
197
198 RecorderOfElectronProcess[CounterOfElectro
199 CounterOfElectronProcess++;
200 PrintInfoDefinition();
201
202 } else {
203
204 RecorderOfPositronProcess[CounterOfPositro
205 CounterOfPositronProcess++;
206 }
207
208 // Build mean free path data using cut value
209
210 if( theMeanFreePath ) delete theMeanFreePath
211 theMeanFreePath = crossSectionHandler->
212 BuildMeanFreePathForMateri
213
214 if(verboseLevel > 0) {
215 G4cout << "The MeanFreePath table is built
216 << G4endl;
217 if(verboseLevel > 1) theMeanFreePath->Prin
218 }
219
220 // Build common DEDX table for all ionisatio
221
222 BuildDEDXTable(aParticleType);
223
224 if (verboseLevel > 0) {
225 G4cout << "G4LowEnergyIonisation::BuildPhy
226 << G4endl;
227 }
228 }
229
230
231 void G4LowEnergyIonisation::BuildLossTable(con
232 {
233 // Build table for energy loss due to soft b
234 // the tables are built for *MATERIALS* binn
235
236 G4double lowKineticEnergy = GetLowerBoundEl
237 G4double highKineticEnergy = GetUpperBoundEl
238 size_t totBin = GetNbinEloss();
239
240 // create table
241
242 if (theLossTable) {
243 theLossTable->clearAndDestroy();
244 delete theLossTable;
245 }
246 const G4ProductionCutsTable* theCoupleTable=
247 G4ProductionCutsTable::GetProductionCu
248 size_t numOfCouples = theCoupleTable->GetTab
249 theLossTable = new G4PhysicsTable(numOfCoupl
250
251 if (shellVacancy != 0) delete shellVacancy;
252 shellVacancy = new G4RDShellVacancy();
253 G4DataVector* ksi = 0;
254 G4DataVector* energy = 0;
255 size_t binForFluo = totBin/10;
256
257 G4PhysicsLogVector* bVector = new G4PhysicsL
258 highKineticEner
259 binForFluo);
260 const G4RDAtomicTransitionManager* transitio
261
262 // Clean up the vector of cuts
263
264 cutForDelta.clear();
265
266 // Loop for materials
267
268 for (size_t m=0; m<numOfCouples; m++) {
269
270 // create physics vector and fill it
271 G4PhysicsLogVector* aVector = new G4Physic
272 highKineticEnergy,
273 totBin);
274
275 // get material parameters needed for the
276 const G4MaterialCutsCouple* couple = theCo
277 const G4Material* material= couple->GetMat
278
279 // the cut cannot be below lowest limit
280 G4double tCut = (*(theCoupleTable->GetEner
281 if(tCut > highKineticEnergy) tCut = highKi
282 cutForDelta.push_back(tCut);
283 const G4ElementVector* theElementVector =
284 size_t NumberOfElements = material->GetNum
285 const G4double* theAtomicNumDensityVector
286 material->GetAtomicNumDens
287 if(verboseLevel > 0) {
288 G4cout << "Energy loss for material # "
289 << " tCut(keV)= " << tCut/keV
290 << G4endl;
291 }
292
293 // now comes the loop for the kinetic ener
294 for (size_t i = 0; i<totBin; i++) {
295
296 G4double lowEdgeEnergy = aVector->GetLow
297 G4double ionloss = 0.;
298
299 // loop for elements in the material
300 for (size_t iel=0; iel<NumberOfElements;
301
302 G4int Z = (G4int)((*theElementVector)[
303
304 G4int nShells = transitionManager->NumberOfS
305
306 for (G4int n=0; n<nShells; n++) {
307
308 G4double e = energySpectrum->Average
309
310 G4double cs= crossSectionHandler->Fi
311 ionloss += e * cs * theAtomicNumDe
312
313 if(verboseLevel > 1 || (Z == 14 && l
314 G4cout << "Z= " << Z
315 << " shell= " << n
316 << " E(keV)= " << lowEdgeEn
317 << " Eav(keV)= " << e/keV
318 << " cs= " << cs
319 << " loss= " << ionloss
320 << " rho= " << theAtomicNumDensit
321 << G4endl;
322 }
323 }
324 G4double esp = energySpectrum->Excitat
325 ionloss += esp * theAtomicNumDensity
326
327 }
328 if(verboseLevel > 1 || (m == 0 && lowEdg
329 G4cout << "Sum: "
330 << " E(keV)= " << lowEdgeEn
331 << " loss(MeV/mm)= " << ionloss*m
332 << G4endl;
333 }
334 aVector->PutValue(i,ionloss);
335 }
336 theLossTable->insert(aVector);
337
338 // fill data for fluorescence
339
340 G4RDVDataSetAlgorithm* interp = new G4RDLo
341 G4RDVEMDataSet* xsis = new G4RDCompositeEM
342 for (size_t iel=0; iel<NumberOfElements; i
343
344 G4int Z = (G4int)((*theElementVector)[ie
345 energy = new G4DataVector();
346 ksi = new G4DataVector();
347
348 for (size_t j = 0; j<binForFluo; j++) {
349
350 G4double lowEdgeEnergy = bVector->GetL
351 G4double cross = 0.;
352 G4double eAverage= 0.;
353 G4int nShells = transitionManager->NumberOfS
354
355 for (G4int n=0; n<nShells; n++) {
356
357 G4double e = energySpectrum->Average
358
359 G4double pro = energySpectrum->Proba
360
361 G4double cs= crossSectionHandler->Fi
362 eAverage += e * cs * theAtomicNumD
363 cross += cs * pro * theAtomicNu
364 if(verboseLevel > 1) {
365 G4cout << "Z= " << Z
366 << " shell= " << n
367 << " E(keV)= " << lowEdgeEn
368 << " Eav(keV)= " << e/keV
369 << " pro= " << pro
370 << " cs= " << cs
371 << G4endl;
372 }
373 }
374
375 G4double coeff = 0.0;
376 if(eAverage > 0.) {
377 coeff = cross/eAverage;
378 eAverage /= cross;
379 }
380
381 if(verboseLevel > 1) {
382 G4cout << "Ksi Coefficient for Z=
383 << " E(keV)= " << lowEdgeEn
384 << " Eav(keV)= " << eAverag
385 << " coeff= " << coeff
386 << G4endl;
387 }
388
389 energy->push_back(lowEdgeEnergy);
390 ksi->push_back(coeff);
391 }
392 interp = new G4RDLogLogInterpolation();
393 G4RDVEMDataSet* set = new G4RDEMDataSet(
394 xsis->AddComponent(set);
395 }
396 if(verboseLevel) xsis->PrintData();
397 shellVacancy->AddXsiTable(xsis);
398 }
399 delete bVector;
400 }
401
402
403 G4VParticleChange* G4LowEnergyIonisation::Post
404 const G4Step& step)
405 {
406 // Delta electron production mechanism on ba
407 // J. Stepanek " A program to determine the
408 // to a single atomic subshell ionisation by
409 // deexcitation or decay of radionuclides",
410 // Comp. Phys. Comm. 1206 pp 1-19 (1997)
411
412 aParticleChange.Initialize(track);
413
414 const G4MaterialCutsCouple* couple = track.G
415 G4double kineticEnergy = track.GetKineticEne
416
417 // Select atom and shell
418
419 G4int Z = crossSectionHandler->SelectRandomA
420 G4int shell = crossSectionHandler->SelectRan
421 const G4RDAtomicShell* atomicShell =
422 (G4RDAtomicTransitionManager::
423 G4double bindingEnergy = atomicShell->Bindin
424 G4int shellId = atomicShell->ShellId();
425
426 // Sample delta energy
427
428 G4int index = couple->GetIndex();
429 G4double tCut = cutForDelta[index];
430 G4double tmax = energySpectrum->MaxEnergyO
431 G4double tDelta = energySpectrum->SampleEner
432
433
434 if(tDelta == 0.0)
435 return G4VContinuousDiscreteProcess::PostS
436
437 // Transform to shell potential
438 G4double deltaKinE = tDelta + 2.0*bindingEne
439 G4double primaryKinE = kineticEnergy + 2.0*b
440
441 // sampling of scattering angle neglecting a
442 G4double deltaMom = std::sqrt(deltaKinE*(del
443 G4double primaryMom = std::sqrt(primaryKinE*
444
445 G4double cost = deltaKinE * (primaryKinE + 2
446 / (deltaMom * prim
447
448 if (cost > 1.) cost = 1.;
449 G4double sint = std::sqrt(1. - cost*cost);
450 G4double phi = twopi * G4UniformRand();
451 G4double dirx = sint * std::cos(phi);
452 G4double diry = sint * std::sin(phi);
453 G4double dirz = cost;
454
455 // Rotate to incident electron direction
456 G4ThreeVector primaryDirection = track.GetMo
457 G4ThreeVector deltaDir(dirx,diry,dirz);
458 deltaDir.rotateUz(primaryDirection);
459 dirx = deltaDir.x();
460 diry = deltaDir.y();
461 dirz = deltaDir.z();
462
463
464 // Take into account atomic motion del is re
465 // kinetic energy of the motion == bindingEn
466
467 cost = 2.0*G4UniformRand() - 1.0;
468 sint = std::sqrt(1. - cost*cost);
469 phi = twopi * G4UniformRand();
470 G4double del = std::sqrt(bindingEnergy *(bin
471 / deltaMom;
472 dirx += del* sint * std::cos(phi);
473 diry += del* sint * std::sin(phi);
474 dirz += del* cost;
475
476 // Find out new primary electron direction
477 G4double finalPx = primaryMom*primaryDirecti
478 G4double finalPy = primaryMom*primaryDirecti
479 G4double finalPz = primaryMom*primaryDirecti
480
481 // create G4DynamicParticle object for delta
482 G4DynamicParticle* theDeltaRay = new G4Dynam
483 theDeltaRay->SetKineticEnergy(tDelta);
484 G4double norm = 1.0/std::sqrt(dirx*dirx + di
485 dirx *= norm;
486 diry *= norm;
487 dirz *= norm;
488 theDeltaRay->SetMomentumDirection(dirx, diry
489 theDeltaRay->SetDefinition(G4Electron::Elect
490
491 G4double theEnergyDeposit = bindingEnergy;
492
493 // fill ParticleChange
494 // changed energy and momentum of the actual
495
496 G4double finalKinEnergy = kineticEnergy - tD
497 if(finalKinEnergy < 0.0) {
498 theEnergyDeposit += finalKinEnergy;
499 finalKinEnergy = 0.0;
500 aParticleChange.ProposeTrackStatus(fStopAn
501
502 } else {
503
504 G4double norm = 1.0/std::sqrt(finalPx*fina
505 finalPx *= norm;
506 finalPy *= norm;
507 finalPz *= norm;
508 aParticleChange.ProposeMomentumDirection(f
509 }
510
511 aParticleChange.ProposeEnergy(finalKinEnergy
512
513 // Generation of Fluorescence and Auger
514 size_t nSecondaries = 0;
515 size_t totalNumber = 1;
516 std::vector<G4DynamicParticle*>* secondaryVe
517 G4DynamicParticle* aSecondary = 0;
518 G4ParticleDefinition* type = 0;
519
520 // Fluorescence data start from element 6
521
522 if (Fluorescence() && Z > 5 && (bindingEnerg
523 || bindingEnergy >= cutForElectro
524
525 secondaryVector = deexcitationManager.Gene
526
527 if (secondaryVector != 0) {
528
529 nSecondaries = secondaryVector->size();
530 for (size_t i = 0; i<nSecondaries; i++)
531
532 aSecondary = (*secondaryVector)[i];
533 if (aSecondary) {
534
535 G4double e = aSecondary->GetKineticE
536 type = aSecondary->GetDefinition();
537 if (e < theEnergyDeposit &&
538 ((type == G4Gamma::Gamma() &&
539 (type == G4Electron::Electron
540
541 theEnergyDeposit -= e;
542 totalNumber++;
543
544 } else {
545
546 delete aSecondary;
547 (*secondaryVector)[i] = 0;
548 }
549 }
550 }
551 }
552 }
553
554 // Save delta-electrons
555
556 aParticleChange.SetNumberOfSecondaries(total
557 aParticleChange.AddSecondary(theDeltaRay);
558
559 // Save Fluorescence and Auger
560
561 if (secondaryVector) {
562
563 for (size_t l = 0; l < nSecondaries; l++)
564
565 aSecondary = (*secondaryVector)[l];
566
567 if(aSecondary) {
568
569 aParticleChange.AddSecondary(aSecondar
570 }
571 }
572 delete secondaryVector;
573 }
574
575 if(theEnergyDeposit < 0.) {
576 G4cout << "G4LowEnergyIonisation: Negative
577 << theEnergyDeposit/eV << " eV" <<
578 theEnergyDeposit = 0.0;
579 }
580 aParticleChange.ProposeLocalEnergyDeposit(th
581
582 return G4VContinuousDiscreteProcess::PostSte
583 }
584
585
586 void G4LowEnergyIonisation::PrintInfoDefinitio
587 {
588 G4String comments = "Total cross sections fr
589 comments += "\n Gamma energy sampled fr
590 comments += "\n Implementation of the c
591 comments += "\n At present it can be us
592 comments += "in the energy range [250eV,100G
593 comments += "\n The process must work w
594
595 G4cout << G4endl << GetProcessName() << ":
596 }
597
598 G4bool G4LowEnergyIonisation::IsApplicable(con
599 {
600 return ( (&particle == G4Electron::Electron
601 }
602
603 std::vector<G4DynamicParticle*>*
604 G4LowEnergyIonisation::DeexciteAtom(const G4Ma
605 G4double incidentEnerg
606 G4double eLoss)
607 {
608 // create vector of secondary particles
609 const G4Material* material = couple->GetMate
610
611 std::vector<G4DynamicParticle*>* partVector
612 new std::vect
613
614 if(eLoss > cutForPhotons && eLoss > cutForEl
615
616 const G4RDAtomicTransitionManager* transit
617 G4RDAtomicTrans
618
619 size_t nElements = material->GetNumberOfEl
620 const G4ElementVector* theElementVector =
621
622 std::vector<G4DynamicParticle*>* secVector
623 G4DynamicParticle* aSecondary = 0;
624 G4ParticleDefinition* type = 0;
625 G4double e;
626 G4ThreeVector position;
627 G4int shell, shellId;
628
629 // sample secondaries
630
631 G4double eTot = 0.0;
632 std::vector<G4int> n =
633 shellVacancy->GenerateNumberOfIonis
634
635 for (size_t i=0; i<nElements; i++) {
636
637 G4int Z = (G4int)((*theElementVector)[i]
638 size_t nVacancies = n[i];
639
640 G4double maxE = transitionManager->Shell
641
642 if (nVacancies && Z > 5 && (maxE>cutForP
643
644 for (size_t j=0; j<nVacancies; j++) {
645
646 shell = crossSectionHandler->SelectRandomS
647 shellId = transitionManager->Shell(Z
648 G4double maxEShell =
649 transitionManager->Shell(
650
651 if (maxEShell>cutForPhotons || maxES
652
653 secVector = deexcitationManager.Generate
654
655 if (secVector != 0) {
656
657 for (size_t l = 0; l<secVector->size()
658
659 aSecondary = (*secVector)[l];
660 if (aSecondary != 0) {
661
662 e = aSecondary->GetKineticEnergy()
663 type = aSecondary->GetDefinition()
664 if ( eTot + e <= eLoss &&
665 ((type == G4Gamma::Gamma() && e
666 (type == G4Electron::Electron()
667
668 eTot += e;
669 partVector->push_bac
670
671 } else {
672
673 delete aSecondary;
674
675 }
676 }
677 }
678 delete secVector;
679 }
680 }
681 }
682 }
683 }
684 }
685 return partVector;
686 }
687
688 G4double G4LowEnergyIonisation::GetMeanFreePat
689 G4double , // previousStepSize
690 G4ForceCondition* cond)
691 {
692 *cond = NotForced;
693 G4int index = (track.GetMaterialCutsCouple(
694 const G4RDVEMDataSet* data = theMeanFreePat
695 G4double meanFreePath = data->FindValue(tra
696 return meanFreePath;
697 }
698
699 void G4LowEnergyIonisation::SetCutForLowEnSecP
700 {
701 cutForPhotons = cut;
702 deexcitationManager.SetCutForSecondaryPhoton
703 }
704
705 void G4LowEnergyIonisation::SetCutForLowEnSecE
706 {
707 cutForElectrons = cut;
708 deexcitationManager.SetCutForAugerElectrons(
709 }
710
711 void G4LowEnergyIonisation::ActivateAuger(G4bo
712 {
713 deexcitationManager.ActivateAugerElectronPro
714 }
715
716