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