source: trunk/source/processes/electromagnetic/lowenergy/src/G4AtomicDeexcitation.cc @ 1350

Last change on this file since 1350 was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4

File size: 15.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//
27// $Id: G4AtomicDeexcitation.cc,v 1.11
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// Authors: Elena Guardincerri (Elena.Guardincerri@ge.infn.it)
31//          Alfonso Mantero (Alfonso.Mantero@ge.infn.it)
32//
33// History:
34// -----------
35// 
36//  16 Sept 2001  First committed to cvs
37//  12 Sep  2003  Bug in auger production fixed
38//
39// -------------------------------------------------------------------
40
41#include "G4AtomicDeexcitation.hh"
42#include "Randomize.hh"
43#include "G4Gamma.hh"
44#include "G4Electron.hh"
45#include "G4AtomicTransitionManager.hh"
46#include "G4FluoTransition.hh"
47
48G4AtomicDeexcitation::G4AtomicDeexcitation():
49  minGammaEnergy(100.*eV),
50  minElectronEnergy(100.*eV),
51  fAuger(false)
52{}
53
54G4AtomicDeexcitation::~G4AtomicDeexcitation()
55{}
56
57std::vector<G4DynamicParticle*>* G4AtomicDeexcitation::GenerateParticles(G4int Z,G4int givenShellId)
58{ 
59
60  std::vector<G4DynamicParticle*>* vectorOfParticles;
61  vectorOfParticles = new std::vector<G4DynamicParticle*>;
62
63  G4DynamicParticle* aParticle;
64  G4int provShellId = 0;
65  G4int counter = 0;
66 
67  // The aim of this loop is to generate more than one fluorecence photon
68  // from the same ionizing event
69  do
70    {
71      if (counter == 0) 
72        // First call to GenerateParticles(...):
73        // givenShellId is given by the process
74        {
75          provShellId = SelectTypeOfTransition(Z, givenShellId);
76         
77          if  ( provShellId >0) 
78            {
79              aParticle = GenerateFluorescence(Z,givenShellId,provShellId); 
80            }
81          else if ( provShellId == -1)
82            {
83              aParticle = GenerateAuger(Z, givenShellId);
84            }
85          else
86            {
87              G4Exception("G4AtomicDeexcitation: starting shell uncorrect: check it");
88            }
89        }
90      else 
91        // Following calls to GenerateParticles(...):
92        // newShellId is given by GenerateFluorescence(...)
93        {
94          provShellId = SelectTypeOfTransition(Z,newShellId);
95          if  (provShellId >0)
96            {
97              aParticle = GenerateFluorescence(Z,newShellId,provShellId);
98            }
99          else if ( provShellId == -1)
100            {
101              aParticle = GenerateAuger(Z, newShellId);
102            }
103          else
104            {
105              G4Exception("G4AtomicDeexcitation: starting shell uncorrect: check it");
106            }
107        }
108      counter++;
109      if (aParticle != 0) {vectorOfParticles->push_back(aParticle);}
110      else {provShellId = -2;}
111    }
112 
113  // Look this in a particular way: only one auger emitted! // ????
114  while (provShellId > -2); 
115
116  // debug 
117  // if (vectorOfParticles->size() > 0) {
118  //   G4cout << " DEEXCITATION!" << G4endl;
119  // }
120
121  return vectorOfParticles;
122}
123
124G4int G4AtomicDeexcitation::SelectTypeOfTransition(G4int Z, G4int shellId)
125{
126  if (shellId <=0 ) 
127    {G4Exception("G4AtomicDeexcitation: zero or negative shellId");}
128
129  G4bool fluoTransitionFoundFlag = false;
130 
131  const G4AtomicTransitionManager*  transitionManager = 
132        G4AtomicTransitionManager::Instance();
133  G4int provShellId = -1;
134  G4int shellNum = 0;
135  G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z); 
136 
137  const G4FluoTransition* refShell = transitionManager->ReachableShell(Z,maxNumOfShells-1);
138
139  // This loop gives shellNum the value of the index of shellId
140  // in the vector storing the list of the shells reachable through
141  // a radiative transition
142  if ( shellId <= refShell->FinalShellId())
143    {
144      while (shellId != transitionManager->ReachableShell(Z,shellNum)->FinalShellId())
145        {
146          if(shellNum ==maxNumOfShells-1)
147            {
148              break;
149            }
150          shellNum++;
151        }
152      G4int transProb = 0; //AM change 29/6/07 was 1
153   
154      G4double partialProb = G4UniformRand();     
155      G4double partSum = 0;
156      const G4FluoTransition* aShell = transitionManager->ReachableShell(Z,shellNum);     
157      G4int trSize =  (aShell->TransitionProbabilities()).size();
158   
159      // Loop over the shells wich can provide an electron for a
160      // radiative transition towards shellId:
161      // in every loop the partial sum of the first transProb shells
162      // is calculated and compared with a random number [0,1].
163      // If the partial sum is greater, the shell whose index is transProb
164      // is chosen as the starting shell for a radiative transition
165      // and its identity is returned
166      // Else, terminateded the loop, -1 is returned
167      while(transProb < trSize){
168       
169         partSum += aShell->TransitionProbability(transProb);
170
171         if(partialProb <= partSum)
172           {
173             provShellId = aShell->OriginatingShellId(transProb);
174             fluoTransitionFoundFlag = true;
175
176             break;
177           }
178         transProb++;
179      }
180
181      // here provShellId is the right one or is -1.
182      // if -1, the control is passed to the Auger generation part of the package
183    }
184
185
186
187  else 
188    {
189 
190     provShellId = -1;
191
192    }
193  return provShellId;
194}
195
196G4DynamicParticle* G4AtomicDeexcitation::GenerateFluorescence(G4int Z, 
197                                                              G4int shellId,
198                                                              G4int provShellId )
199{ 
200
201
202  const G4AtomicTransitionManager*  transitionManager = G4AtomicTransitionManager::Instance();
203  //  G4int provenienceShell = provShellId;
204
205  //isotropic angular distribution for the outcoming photon
206  G4double newcosTh = 1.-2.*G4UniformRand();
207  G4double  newsinTh = std::sqrt(1.-newcosTh*newcosTh);
208  G4double newPhi = twopi*G4UniformRand();
209 
210  G4double xDir =  newsinTh*std::sin(newPhi);
211  G4double yDir = newsinTh*std::cos(newPhi);
212  G4double zDir = newcosTh;
213 
214  G4ThreeVector newGammaDirection(xDir,yDir,zDir);
215 
216  G4int shellNum = 0;
217  G4int maxNumOfShells = transitionManager->NumberOfReachableShells(Z);
218 
219  // find the index of the shell named shellId
220  while (shellId != transitionManager->
221         ReachableShell(Z,shellNum)->FinalShellId())
222    {
223      if(shellNum == maxNumOfShells-1)
224        {
225          break;
226        }
227      shellNum++;
228    }
229  // number of shell from wich an electron can reach shellId
230  size_t transitionSize = transitionManager->
231    ReachableShell(Z,shellNum)->OriginatingShellIds().size();
232 
233  size_t index = 0;
234 
235  // find the index of the shell named provShellId in the vector
236  // storing the shells from which shellId can be reached
237  while (provShellId != transitionManager->
238         ReachableShell(Z,shellNum)->OriginatingShellId(index))
239    {
240      if(index ==  transitionSize-1)
241        {
242          break;
243        }
244      index++;
245    }
246  // energy of the gamma leaving provShellId for shellId
247  G4double transitionEnergy = transitionManager->
248    ReachableShell(Z,shellNum)->TransitionEnergy(index);
249 
250  // This is the shell where the new vacancy is: it is the same
251  // shell where the electron came from
252  newShellId = transitionManager->
253    ReachableShell(Z,shellNum)->OriginatingShellId(index);
254 
255 
256  G4DynamicParticle* newPart = new G4DynamicParticle(G4Gamma::Gamma(), 
257                                                     newGammaDirection,
258                                                     transitionEnergy);
259  return newPart;
260}
261
262G4DynamicParticle* G4AtomicDeexcitation::GenerateAuger(G4int Z, G4int shellId)
263{
264  if(!fAuger) return 0;
265 
266
267  const G4AtomicTransitionManager*  transitionManager = 
268        G4AtomicTransitionManager::Instance();
269
270
271
272  if (shellId <=0 ) 
273    {G4Exception("G4AtomicDeexcitation: zero or negative shellId");}
274 
275  // G4int provShellId = -1;
276  G4int maxNumOfShells = transitionManager->NumberOfReachableAugerShells(Z); 
277 
278  const G4AugerTransition* refAugerTransition = 
279        transitionManager->ReachableAugerShell(Z,maxNumOfShells-1);
280
281
282  // This loop gives to shellNum the value of the index of shellId
283  // in the vector storing the list of the vacancies in the variuos shells
284  // that can originate a NON-radiative transition
285 
286  // ---- MGP ---- Next line commented out to remove compilation warning
287  // G4int p = refAugerTransition->FinalShellId();
288
289  G4int shellNum = 0;
290
291
292  if ( shellId <= refAugerTransition->FinalShellId() ) 
293    //"FinalShellId" is final from the point of view of the elctron who makes the transition,
294    // being the Id of the shell in which there is a vacancy
295    {
296      G4int pippo = transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId();
297      if (shellId  != pippo ) {
298        do { 
299          shellNum++;
300          if(shellNum == maxNumOfShells)
301            {
302
303              //G4Exception("G4AtomicDeexcitation: No Auger transition found");
304              return 0;
305            }
306        }
307        while (shellId != (transitionManager->ReachableAugerShell(Z,shellNum)->FinalShellId()) ) ;
308      }
309
310
311      // Now we have that shellnum is the shellIndex of the shell named ShellId
312
313      //      G4cout << " the index of the shell is: "<<shellNum<<G4endl;
314
315      // But we have now to select two shells: one for the transition,
316      // and another for the auger emission.
317
318      G4int transitionLoopShellIndex = 0;     
319      G4double partSum = 0;
320      const G4AugerTransition* anAugerTransition = 
321            transitionManager->ReachableAugerShell(Z,shellNum);
322
323      //      G4cout << " corresponding to the ID: "<< anAugerTransition->FinalShellId() << G4endl;
324
325
326      G4int transitionSize = 
327            (anAugerTransition->TransitionOriginatingShellIds())->size();
328      while (transitionLoopShellIndex < transitionSize) {
329
330        std::vector<G4int>::const_iterator pos = 
331               anAugerTransition->TransitionOriginatingShellIds()->begin();
332
333        G4int transitionLoopShellId = *(pos+transitionLoopShellIndex);
334        G4int numberOfPossibleAuger = 
335              (anAugerTransition->AugerTransitionProbabilities(transitionLoopShellId))->size();
336        G4int augerIndex = 0;
337        //      G4int partSum2 = 0;
338
339
340        if (augerIndex < numberOfPossibleAuger) {
341         
342          do 
343            {
344              G4double thisProb = anAugerTransition->AugerTransitionProbability(augerIndex, 
345                                                                                transitionLoopShellId);
346              partSum += thisProb;
347              augerIndex++;
348             
349            } while (augerIndex < numberOfPossibleAuger);
350                }
351        transitionLoopShellIndex++;
352      }
353     
354
355
356      // Now we have the entire probability of an auger transition for the vacancy
357      // located in shellNum (index of shellId)
358
359      // AM *********************** F I X E D **************************** AM
360      // Here we duplicate the previous loop, this time looking to the sum of the probabilities
361      // to be under the random number shoot by G4 UniformRdandom. This could have been done in the
362      // previuos loop, while integrating the probabilities. There is a bug that will be fixed
363      // 5 minutes from now: a line:
364      // G4int numberOfPossibleAuger = (anAugerTransition->
365      // AugerTransitionProbabilities(transitionLoopShellId))->size();
366      // to be inserted.
367      // AM *********************** F I X E D **************************** AM
368
369      // Remains to get the same result with a single loop.
370
371      // AM *********************** F I X E D **************************** AM
372      // Another Bug: in EADL Auger Transition are normalized to all the transitions deriving from
373      // a vacancy in one shell, but not all of these are present in data tables. So if a transition
374      // doesn't occur in the main one a local energy deposition must occur, instead of (like now)
375      // generating the last transition present in EADL data.
376      // AM *********************** F I X E D **************************** AM
377
378
379      G4double totalVacancyAugerProbability = partSum;
380
381
382      //And now we start to select the right auger transition and emission
383      G4int transitionRandomShellIndex = 0;
384      G4int transitionRandomShellId = 1;
385      G4int augerIndex = 0;
386      partSum = 0; 
387      G4double partialProb = G4UniformRand();
388      // G4int augerOriginatingShellId = 0;
389     
390      G4int numberOfPossibleAuger = 0;
391     
392      G4bool foundFlag = false;
393
394      while (transitionRandomShellIndex < transitionSize) {
395
396        std::vector<G4int>::const_iterator pos = 
397               anAugerTransition->TransitionOriginatingShellIds()->begin();
398
399        transitionRandomShellId = *(pos+transitionRandomShellIndex);
400       
401        augerIndex = 0;
402        numberOfPossibleAuger = (anAugerTransition-> 
403                                 AugerTransitionProbabilities(transitionRandomShellId))->size();
404
405        while (augerIndex < numberOfPossibleAuger) {
406          G4double thisProb =anAugerTransition->AugerTransitionProbability(augerIndex, 
407                                                                           transitionRandomShellId);
408
409          partSum += thisProb;
410         
411          if (partSum >= (partialProb*totalVacancyAugerProbability) ) { // was /
412            foundFlag = true;
413            break;
414          }
415          augerIndex++;
416        }
417        if (partSum >= (partialProb*totalVacancyAugerProbability) ) {break;} // was /
418        transitionRandomShellIndex++;
419      }
420
421      // Now we have the index of the shell from wich comes the auger electron (augerIndex),
422      // and the id of the shell, from which the transition e- come (transitionRandomShellid)
423      // If no Transition has been found, 0 is returned. 
424
425      if (!foundFlag) {return 0;}     
426     
427      // Isotropic angular distribution for the outcoming e-
428      G4double newcosTh = 1.-2.*G4UniformRand();
429      G4double  newsinTh = std::sqrt(1.-newcosTh*newcosTh);
430      G4double newPhi = twopi*G4UniformRand();
431     
432      G4double xDir =  newsinTh*std::sin(newPhi);
433      G4double yDir = newsinTh*std::cos(newPhi);
434      G4double zDir = newcosTh;
435     
436      G4ThreeVector newElectronDirection(xDir,yDir,zDir);
437     
438      // energy of the auger electron emitted
439     
440     
441      G4double transitionEnergy = anAugerTransition->AugerTransitionEnergy(augerIndex, transitionRandomShellId);
442      /*
443        G4cout << "AUger TransitionId " << anAugerTransition->FinalShellId() << G4endl;
444        G4cout << "augerIndex: " << augerIndex << G4endl;
445        G4cout << "transitionShellId: " << transitionRandomShellId << G4endl;
446      */
447     
448      // This is the shell where the new vacancy is: it is the same
449      // shell where the electron came from
450      newShellId = transitionRandomShellId;
451     
452     
453      G4DynamicParticle* newPart = new G4DynamicParticle(G4Electron::Electron(), 
454                                                         newElectronDirection,
455                                                         transitionEnergy);
456      return newPart;
457
458    }
459  else 
460    {
461      //G4Exception("G4AtomicDeexcitation: no auger transition found");
462      return 0;
463    }
464 
465}
466
467void G4AtomicDeexcitation::SetCutForSecondaryPhotons(G4double cut)
468{
469  minGammaEnergy = cut;
470}
471
472void G4AtomicDeexcitation::SetCutForAugerElectrons(G4double cut)
473{
474  minElectronEnergy = cut;
475}
476
477void G4AtomicDeexcitation::ActivateAugerElectronProduction(G4bool val)
478{
479  fAuger = val;
480}
481
482
483
484
485
486
487
Note: See TracBrowser for help on using the repository browser.