| 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 // Geant4 Class file
30 //
31 // Authors: Alfonso Mantero (Alfonso.Mantero@g
32 //
33 // Created 22 April 2010 from old G4UAtomicDee
34 //
35 // Modified:
36 // ---------
37 // 20 Oct 2011 Alf modified to take into acc
38 // 03 Nov 2011 Alf Extended Empirical and Fo
39 // out thei ranges with Anal
40 // 07 Nov 2011 Alf Restored original ioniati
41 // letting scaled ones for o
42 // 20 Mar 2012 LP Register G4PenelopeIonisa
43 //
44 // -------------------------------------------
45 //
46 // Class description:
47 // Implementation of atomic deexcitation
48 //
49 // -------------------------------------------
50
51 #include "G4UAtomicDeexcitation.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "Randomize.hh"
55 #include "G4Gamma.hh"
56 #include "G4AtomicTransitionManager.hh"
57 #include "G4FluoTransition.hh"
58 #include "G4Electron.hh"
59 #include "G4Positron.hh"
60 #include "G4Proton.hh"
61 #include "G4Alpha.hh"
62
63 #include "G4teoCrossSection.hh"
64 #include "G4empCrossSection.hh"
65 #include "G4PenelopeIonisationCrossSection.hh"
66 #include "G4LivermoreIonisationCrossSection.hh
67 #include "G4EmCorrections.hh"
68 #include "G4LossTableManager.hh"
69 #include "G4EmParameters.hh"
70 #include "G4Material.hh"
71 #include "G4AtomicShells.hh"
72
73 using namespace std;
74
75 //....oooOO0OOooo........oooOO0OOooo........oo
76
77 G4UAtomicDeexcitation::G4UAtomicDeexcitation()
78 G4VAtomDeexcitation("UAtomDeexcitation"),
79 minGammaEnergy(DBL_MAX),
80 minElectronEnergy(DBL_MAX),
81 newShellId(-1)
82 {
83 anaPIXEshellCS = nullptr;
84 PIXEshellCS = nullptr;
85 ePIXEshellCS = nullptr;
86 emcorr = G4LossTableManager::Instance()->EmC
87 theElectron = G4Electron::Electron();
88 thePositron = G4Positron::Positron();
89 transitionManager = G4AtomicTransitionManage
90 }
91
92 //....oooOO0OOooo........oooOO0OOooo........oo
93
94 G4UAtomicDeexcitation::~G4UAtomicDeexcitation(
95 {
96 delete anaPIXEshellCS;
97 delete PIXEshellCS;
98 delete ePIXEshellCS;
99 }
100
101 //....oooOO0OOooo........oooOO0OOooo........oo
102
103 void G4UAtomicDeexcitation::InitialiseForNewRu
104 {
105 if(!IsFluoActive()) { return; }
106 transitionManager->Initialise();
107 if(!IsPIXEActive()) { return; }
108
109 if(!anaPIXEshellCS) {
110 anaPIXEshellCS = new G4teoCrossSection("EC
111 }
112 G4cout << G4endl;
113 G4cout << "### === G4UAtomicDeexcitation::In
114
115 G4EmParameters* param = G4EmParameters::Inst
116 G4String namePIXExsModel = param->PIXECrossS
117 G4String namePIXExsElectronModel = param->PI
118
119 // Check if old cross section for p/ion shou
120 if(PIXEshellCS && namePIXExsModel != PIXEshe
121 {
122 delete PIXEshellCS;
123 PIXEshellCS = nullptr;
124 }
125
126 // Instantiate new proton/ion cross section
127 if(!PIXEshellCS) {
128 if (namePIXExsModel == "ECPSSR_FormFactor"
129 {
130 PIXEshellCS = new G4teoCrossSection(namePIXE
131 }
132 else if(namePIXExsModel == "ECPSSR_ANSTO")
133 {
134 PIXEshellCS = new G4teoCrossSection(namePIXE
135 }
136 else if(namePIXExsModel == "Empirical")
137 {
138 PIXEshellCS = new G4empCrossSection(namePIXE
139 }
140 }
141
142 // Check if old cross section for e+- should
143 if(ePIXEshellCS && namePIXExsElectronModel !
144 {
145 delete ePIXEshellCS;
146 ePIXEshellCS = nullptr;
147 }
148
149 // Instantiate new e+- cross section
150 if(nullptr == ePIXEshellCS)
151 {
152 if(namePIXExsElectronModel == "Empirical
153 {
154 ePIXEshellCS = new G4empCrossSection("Empi
155 }
156 else if(namePIXExsElectronModel == "ECPS
157 {
158 ePIXEshellCS = new G4teoCrossSection("ECPS
159 }
160 else if (namePIXExsElectronModel == "Pen
161 {
162 ePIXEshellCS = new G4PenelopeIonisationCro
163 }
164 else
165 {
166 ePIXEshellCS = new G4LivermoreIonisationCr
167 }
168 }
169 }
170
171 //....oooOO0OOooo........oooOO0OOooo........oo
172
173 void G4UAtomicDeexcitation::InitialiseForExtra
174 {}
175
176 //....oooOO0OOooo........oooOO0OOooo........oo
177
178 const G4AtomicShell*
179 G4UAtomicDeexcitation::GetAtomicShell(G4int Z,
180 {
181 return transitionManager->Shell(Z, (std::siz
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oo
185
186 void G4UAtomicDeexcitation::GenerateParticles(
187 std::vector<G4DynamicParticle*
188 const G4AtomicShell* atomicShell,
189 G4int Z,
190 G4double gammaCut,
191 G4double eCut)
192 {
193 // Defined initial conditions
194 G4int givenShellId = atomicShell->ShellId();
195 minGammaEnergy = gammaCut;
196 minElectronEnergy = eCut;
197 vacancyArray.clear();
198
199 // generation secondaries
200 G4DynamicParticle* aParticle=0;
201 G4int provShellId = 0;
202
203 //ORIGINAL METHOD BY ALFONSO MANTERO
204 if (!IsAugerCascadeActive())
205 {
206 //----------------------------
207 G4int counter = 0;
208
209 // limits of the EPDL data
210 if (Z>5 && Z<105) {
211
212 // The aim of this loop is to generate more
213 // from the same ionizing event
214 do
215 {
216 if (counter == 0)
217 // First call to GenerateParticles(...
218 // givenShellId is given by the proces
219 {
220 provShellId = SelectTypeOfTransition(Z, gi
221
222 if (provShellId >0)
223 {
224 aParticle =
225 GenerateFluorescence(Z, givenShellId
226 }
227 else if (provShellId == -1)
228 {
229 aParticle = GenerateAuger(Z, givenShel
230 }
231 }
232 else
233 // Following calls to GenerateParticle
234 // newShellId is given by GenerateFluo
235 {
236 provShellId = SelectTypeOfTransition(Z,new
237 if (provShellId >0)
238 {
239 aParticle = GenerateFluorescence(Z,new
240 }
241 else if ( provShellId == -1)
242 {
243 aParticle = GenerateAuger(Z, newShellI
244 }
245 }
246 ++counter;
247 if (aParticle != 0)
248 {
249 vectorOfParticles->push_back(aParticle);
250 }
251 else {provShellId = -2;}
252 }
253 while (provShellId > -2);
254 }
255 } // Auger cascade is not active
256
257 //END OF ORIGINAL METHOD BY ALFONSO MANTERO
258 //----------------------
259
260 // NEW METHOD
261 // Auger cascade by Burkhant Suerfu on March
262 if (IsAugerCascadeActive())
263 {
264 //----------------------
265 vacancyArray.push_back(givenShellId);
266
267 // let's check that 5<Z<100
268 if (Z<6 || Z>104){
269 return;
270 }
271
272 // as long as there is vacancy to be fil
273 while(!vacancyArray.empty()){
274 // prepare to process the last element, and
275 givenShellId = vacancyArray[0];
276 provShellId = SelectTypeOfTransition(Z,given
277
278 //G4cout<<"\n------ Atom Transition with Z:
279 // <<givenShellId<<" & target:"<<provShel
280 if(provShellId>0){
281 aParticle = GenerateFluorescence(Z,givenSh
282 }
283 else if(provShellId == -1){
284 aParticle = GenerateAuger(Z, givenShellId)
285 }
286 // if a particle is created, put it in the
287 if(aParticle!=0)
288 vectorOfParticles->push_back(aParticle);
289
290 // one vacancy has been processed. Erase it
291 vacancyArray.erase(vacancyArray.begin());
292 }
293 //----------------------
294 //End of Auger cascade by Burkhant Suerf
295
296 } // Auger cascade is active
297 }
298
299 //....oooOO0OOooo........oooOO0OOooo........oo
300
301 G4double
302 G4UAtomicDeexcitation::GetShellIonisationCross
303 const G4ParticleDefinition* pdef,
304 G4int Z,
305 G4AtomicShellEnumerator shellEnum,
306 G4double kineticEnergy,
307 const G4Material* mat)
308 {
309 // we must put a control on the shell that a
310 // some shells should not pass (line "0" or
311
312 // check atomic number
313 G4double xsec = 0.0;
314 if(Z > 93 || Z < 6 ) { return xsec; } //corr
315 G4int idx = G4int(shellEnum);
316 if(idx >= G4AtomicShells::GetNumberOfShells(
317
318 if(pdef == theElectron || pdef == thePositro
319 xsec = ePIXEshellCS->CrossSection(Z,shellE
320 return xsec;
321 }
322
323 G4double mass = pdef->GetPDGMass();
324 G4double escaled = kineticEnergy;
325 G4double q2 = 0.0;
326
327 // scaling to protons for all particles excl
328 G4int pdg = pdef->GetPDGEncoding();
329 if (pdg != 2212 && pdg != 1000020040)
330 {
331 mass = proton_mass_c2;
332 escaled = kineticEnergy*mass/(pdef->GetP
333
334 if(mat) {
335 q2 = emcorr->EffectiveChargeSquareRatio(pdef
336 } else {
337 G4double q = pdef->GetPDGCharge()/eplus;
338 q2 = q*q;
339 }
340 }
341
342 if(PIXEshellCS) {
343 xsec = PIXEshellCS->CrossSection(Z,shellEn
344 }
345 if(xsec < 1e-100) {
346 xsec = anaPIXEshellCS->CrossSection(Z,shel
347 }
348
349 if (q2) {xsec *= q2;}
350
351 return xsec;
352 }
353
354 //....oooOO0OOooo........oooOO0OOooo........oo
355
356 void G4UAtomicDeexcitation::SetCutForSecondary
357 {
358 minGammaEnergy = cut;
359 }
360
361 //....oooOO0OOooo........oooOO0OOooo........oo
362
363 void G4UAtomicDeexcitation::SetCutForAugerElec
364 {
365 minElectronEnergy = cut;
366 }
367
368 //....oooOO0OOooo........oooOO0OOooo........oo
369
370 G4double G4UAtomicDeexcitation::ComputeShellIo
371 const G4ParticleDefinition* p,
372 G4int Z,
373 G4AtomicShellEnumerator shell,
374 G4double kinE,
375 const G4Material* mat)
376 {
377 return GetShellIonisationCrossSectionPerAtom
378 }
379
380 //....oooOO0OOooo........oooOO0OOooo........oo
381
382 G4int G4UAtomicDeexcitation::SelectTypeOfTrans
383 {
384 if (shellId <=0 ) {
385 return 0;
386 }
387
388 G4int provShellId = -1;
389 G4int shellNum = 0;
390 G4int maxNumOfShells = transitionManager->Nu
391
392 const G4FluoTransition* refShell =
393 transitionManager->ReachableShell(Z,maxNum
394
395 // This loop gives shellNum the value of the
396 // in the vector storing the list of the she
397 // a radiative transition
398 if ( shellId <= refShell->FinalShellId())
399 {
400 while (shellId != transitionManager->Rea
401 {
402 if(shellNum ==maxNumOfShells-1)
403 {
404 break;
405 }
406 shellNum++;
407 }
408 G4int transProb = 0; //AM change 29/6/07
409
410 G4double partialProb = G4UniformRand();
411 G4double partSum = 0;
412 const G4FluoTransition* aShell = transit
413 G4int trSize = (G4int)(aShell->Transiti
414
415 // Loop over the shells wich can provide
416 // radiative transition towards shellId:
417 // in every loop the partial sum of the
418 // is calculated and compared with a ran
419 // If the partial sum is greater, the sh
420 // is chosen as the starting shell for a
421 // and its identity is returned
422 // Else, terminateded the loop, -1 is re
423 while(transProb < trSize){
424 partSum += aShell->TransitionProbability(tra
425
426 if(partialProb <= partSum)
427 {
428 provShellId = aShell->OriginatingShellId
429 break;
430 }
431 ++transProb;
432 }
433 // here provShellId is the right one or
434 // if -1, the control is passed to the A
435 }
436 else
437 {
438 provShellId = -1;
439 }
440 return provShellId;
441 }
442
443 //....oooOO0OOooo........oooOO0OOooo........oo
444
445 G4DynamicParticle*
446 G4UAtomicDeexcitation::GenerateFluorescence(G4
447 G4int provShellId )
448 {
449 if (shellId <=0 )
450 {
451 return nullptr;
452 }
453
454 //isotropic angular distribution for the out
455 G4double newcosTh = 1.-2.*G4UniformRand();
456 G4double newsinTh = std::sqrt((1.-newcosTh)*
457 G4double newPhi = twopi*G4UniformRand();
458
459 G4double xDir = newsinTh*std::sin(newPhi);
460 G4double yDir = newsinTh*std::cos(newPhi);
461 G4double zDir = newcosTh;
462
463 G4ThreeVector newGammaDirection(xDir,yDir,zD
464
465 G4int shellNum = 0;
466 G4int maxNumOfShells = transitionManager->Nu
467
468 // find the index of the shell named shellId
469 while (shellId != transitionManager->
470 ReachableShell(Z,shellNum)->FinalShellId())
471 {
472 if(shellNum == maxNumOfShells-1)
473 {
474 break;
475 }
476 ++shellNum;
477 }
478 // number of shell from wich an electron can
479 G4int transitionSize = (G4int)transitionMana
480 ReachableShell(Z,shellNum)->OriginatingShe
481
482 G4int index = 0;
483
484 // find the index of the shell named provShe
485 // storing the shells from which shellId can
486 while (provShellId != transitionManager->
487 ReachableShell(Z,shellNum)->OriginatingShel
488 {
489 if(index == transitionSize-1)
490 {
491 break;
492 }
493 ++index;
494 }
495 // energy of the gamma leaving provShellId f
496 G4double transitionEnergy = transitionManage
497 ReachableShell(Z,shellNum)->TransitionEner
498
499 if (transitionEnergy < minGammaEnergy) retur
500
501 // This is the shell where the new vacancy i
502 // shell where the electron came from
503 newShellId = transitionManager->
504 ReachableShell(Z,shellNum)->OriginatingShe
505
506 G4DynamicParticle* newPart = new G4DynamicPa
507 newGammaDirection,
508 transitionEnergy);
509
510 //Auger cascade by Burkhant Suerfu on March
511 if (IsAugerCascadeActive()) vacancyArray.pus
512
513 return newPart;
514 }
515
516 //....oooOO0OOooo........oooOO0OOooo........oo
517
518 G4DynamicParticle* G4UAtomicDeexcitation::Gene
519 {
520 if(!IsAugerActive()) {
521 // G4cout << "auger inactive!" << G4end
522 return nullptr;
523 }
524
525 if (shellId <=0 ) {
526 //G4Exception("G4UAtomicDeexcitation::Gene
527 // JustWarning, "Energy deposited local
528 return nullptr;
529 }
530
531 G4int maxNumOfShells = transitionManager->Nu
532
533 const G4AugerTransition* refAugerTransition
534 transitionManager->ReachableAugerShell(Z,m
535
536 // This loop gives to shellNum the value of
537 // in the vector storing the list of the vac
538 // that can originate a NON-radiative transi
539 G4int shellNum = 0;
540
541 if ( shellId <= refAugerTransition->FinalShe
542 // "FinalShellId" is final from the point
543 // who makes the transition,
544 // being the Id of the shell in which ther
545 {
546 G4int pippo = transitionManager->Reachab
547 if (shellId != pippo ) {
548 do {
549 ++shellNum;
550 if(shellNum == maxNumOfShells)
551 {
552 // G4cout << "No Auger transition foun
553 return 0;
554 }
555 }
556 while (shellId != (transitionManager->Reacha
557 }
558
559 // Now we have that shellnum is the shel
560 // G4cout << " the index of the she
561 // But we have now to select two shells:
562 // and another for the auger emission.
563 G4int transitionLoopShellIndex = 0;
564 G4double partSum = 0;
565 const G4AugerTransition* anAugerTransiti
566 transitionManager->ReachableAugerShell(Z,she
567
568 G4int transitionSize = (G4int)
569 (anAugerTransition->TransitionOriginatingShe
570 while (transitionLoopShellIndex < transi
571
572 std::vector<G4int>::const_iterator pos
573 anAugerTransition->TransitionOriginatingSh
574
575 G4int transitionLoopShellId = *(pos+tr
576 G4int numberOfPossibleAuger = (G4int)
577 (anAugerTransition->AugerTransitionProbabi
578 G4int augerIndex = 0;
579
580 if (augerIndex < numberOfPossibleAuger) {
581 do
582 {
583 G4double thisProb = anAugerTransition-
584 transitionLoopShellId);
585 partSum += thisProb;
586 augerIndex++;
587
588 } while (augerIndex < numberOfPossibleAu
589 }
590 ++transitionLoopShellIndex;
591 }
592
593 G4double totalVacancyAugerProbability =
594
595 //And now we start to select the right a
596 G4int transitionRandomShellIndex = 0;
597 G4int transitionRandomShellId = 1;
598 G4int augerIndex = 0;
599 partSum = 0;
600 G4double partialProb = G4UniformRand();
601
602 G4int numberOfPossibleAuger = 0;
603 G4bool foundFlag = false;
604
605 while (transitionRandomShellIndex < tran
606
607 std::vector<G4int>::const_iterator pos
608 anAugerTransition->TransitionOriginatingSh
609
610 transitionRandomShellId = *(pos+transi
611
612 augerIndex = 0;
613 numberOfPossibleAuger = (G4int)(anAugerTrans
614 AugerTransitionProbabilities(transiti
615
616 while (augerIndex < numberOfPossibleAu
617 G4double thisProb =anAugerTransition->Auge
618 transitionRandomShellId);
619
620 partSum += thisProb;
621
622 if (partSum >= (partialProb*totalVac
623 foundFlag = true;
624 break;
625 }
626 augerIndex++;
627 }
628 if (partSum >= (partialProb*totalVacan
629 ++transitionRandomShellIndex;
630 }
631
632 // Now we have the index of the shell fr
633 // and the id of the shell, from which t
634 // If no Transition has been found, 0 is
635 if (!foundFlag) {
636 return nullptr;
637 }
638
639 // Isotropic angular distribution for th
640 G4double newcosTh = 1.-2.*G4UniformRand(
641 G4double newsinTh = std::sqrt(1.-newcosT
642 G4double newPhi = twopi*G4UniformRand();
643
644 G4double xDir = newsinTh*std::sin(newPhi
645 G4double yDir = newsinTh*std::cos(newPhi
646 G4double zDir = newcosTh;
647
648 G4ThreeVector newElectronDirection(xDir,
649
650 // energy of the auger electron emitted
651 G4double transitionEnergy =
652 anAugerTransition->AugerTransitionEnergy(aug
653
654 if (transitionEnergy < minElectronEnergy
655 return nullptr;
656 }
657
658 // This is the shell where the new vacan
659 // shell where the electron came from
660 newShellId = transitionRandomShellId;
661
662 //Auger cascade by Burkhant Suerfu on Ma
663 if (IsAugerCascadeActive())
664 {
665 vacancyArray.push_back(newShellId);
666 vacancyArray.push_back(anAugerTransition->
667 }
668
669 return new G4DynamicParticle(G4Electron:
670 newElectronDirection,
671 transitionEnergy);
672 }
673 else
674 {
675 return nullptr;
676 }
677 }
678