ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ECDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4ECDecay.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 //
27 // //
28 // File: G4ECDecay.cc //
29 // Author: D.H. Wright (SLAC) //
30 // Date: 25 November 2014 //
31 // //
33 
34 #include "G4ECDecay.hh"
35 #include "G4IonTable.hh"
36 #include "Randomize.hh"
37 #include "G4ThreeVector.hh"
38 #include "G4DynamicParticle.hh"
39 #include "G4DecayProducts.hh"
40 #include "G4VAtomDeexcitation.hh"
41 #include "G4AtomicShells.hh"
42 #include "G4Electron.hh"
43 #include "G4LossTableManager.hh"
44 #include "G4PhysicalConstants.hh"
45 #include "G4SystemOfUnits.hh"
46 
47 G4ECDecay::G4ECDecay(const G4ParticleDefinition* theParentNucleus,
48  const G4double& branch, const G4double& Qvalue,
49  const G4double& excitationE,
50  const G4Ions::G4FloatLevelBase& flb,
51  const G4RadioactiveDecayMode& mode)
52  : G4NuclearDecay("electron capture", mode, excitationE, flb), transitionQ(Qvalue),
53  applyARM(true)
54 {
55  SetParent(theParentNucleus); // Store name of parent nucleus, delete G4MT_parent
56  SetBR(branch);
57 
58  SetNumberOfDaughters(2);
59  G4IonTable* theIonTable =
60  (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
61  G4int daughterZ = theParentNucleus->GetAtomicNumber() - 1;
62  G4int daughterA = theParentNucleus->GetAtomicMass();
63  SetDaughter(0, theIonTable->GetIon(daughterZ, daughterA, excitationE, flb) );
64  SetDaughter(1, "nu_e");
65  DefineSubshellProbabilities(daughterZ, daughterZ);
66 }
67 
68 
69 G4ECDecay::~G4ECDecay()
70 {}
71 
72 
73 G4DecayProducts* G4ECDecay::DecayIt(G4double)
74 {
75  // Fill G4MT_parent with theParentNucleus (stored by SetParent in ctor)
76  CheckAndFillParent();
77 
78  // Fill G4MT_daughters with alpha and residual nucleus (stored by SetDaughter)
79  CheckAndFillDaughters();
80 
81  // Get shell number of captured electron
82  G4int shellIndex = -1;
83  G4double ran;
84  switch (theMode)
85  {
86  case KshellEC:
87  shellIndex = 0;
88  break;
89  case LshellEC: // PL1+PL2+PL3=1
90  ran=G4UniformRand();
91  if (ran <= PL1) shellIndex =1;
92  else if (ran<= (PL1+PL2)) shellIndex =2;
93  else shellIndex =3;
94  break;
95  case MshellEC: // PM1+PM2+PM3=1
96  ran=G4UniformRand();
97  if (ran < PM1) shellIndex =4;
98  else if (ran< (PM1+PM2)) shellIndex =5;
99  else shellIndex = 6;
100  break;
101  case NshellEC: // PN1+PN2+PN3=1
102  ran=G4UniformRand();
103  if (ran < PN1) shellIndex = 9;
104  else if (ran<= (PN1+PN2)) shellIndex =2;
105  else shellIndex =10;
106  break;
107  default:
108  G4Exception("G4ECDecay::DecayIt()", "HAD_RDM_009",
109  FatalException, "Invalid electron shell selected");
110  }
111 
112  // Initialize decay products with parent nucleus at rest
113  G4DynamicParticle parentParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
114  G4DecayProducts* products = new G4DecayProducts(parentParticle);
115  G4double eBind = 0.0;
116 
117  // G4LossTableManager must already be initialized with G4UAtomicDeexcitation
118  // This is currently done in G4RadioactiveDecay::BuildPhysicsTable
119  G4VAtomDeexcitation* atomDeex =
120  G4LossTableManager::Instance()->AtomDeexcitation();
121  std::vector<G4DynamicParticle*> armProducts;
122 
123  if (applyARM) {
124  if (atomDeex) {
125  G4int aZ = G4MT_daughters[0]->GetAtomicNumber();
126  G4int nShells = G4AtomicShells::GetNumberOfShells(aZ);
127  if (shellIndex >= nShells) shellIndex = nShells;
128  G4AtomicShellEnumerator as = G4AtomicShellEnumerator(shellIndex);
129  const G4AtomicShell* shell = atomDeex->GetAtomicShell(aZ, as);
130  eBind = shell->BindingEnergy();
131  if (atomDeex->IsFluoActive() && aZ > 5 && aZ < 100) {
132  // Do atomic relaxation
133  // VI, SI
134  // Allows fixing of Bugzilla 1727
135  //const G4double deexLimit = 0.1*keV;
136  G4double deexLimit = 0.1*keV;
137  if (G4EmParameters::Instance()->DeexcitationIgnoreCut()) deexLimit =0.;
138  //
139  atomDeex->GenerateParticles(&armProducts, shell, aZ, deexLimit, deexLimit);
140  }
141 
142  G4double productEnergy = 0.;
143  for (G4int i = 0; i < G4int(armProducts.size()); i++)
144  productEnergy += armProducts[i]->GetKineticEnergy();
145 
146  G4double deficit = shell->BindingEnergy() - productEnergy;
147  if (deficit > 0.0) {
148  // Add a dummy electron to make up extra energy
149  G4double cosTh = 1.-2.*G4UniformRand();
150  G4double sinTh = std::sqrt(1.- cosTh*cosTh);
151  G4double phi = twopi*G4UniformRand();
152 
153  G4ThreeVector electronDirection(sinTh*std::sin(phi),
154  sinTh*std::cos(phi), cosTh);
155  G4DynamicParticle* extra =
156  new G4DynamicParticle(G4Electron::Electron(), electronDirection,
157  deficit);
158  armProducts.push_back(extra);
159  }
160  } // atomDeex
161  } // applyARM
162 
163  G4double daughterMass = G4MT_daughters[0]->GetPDGMass();
164 
165  // CM momentum using Q value corrected for binding energy of captured electron
166  G4double Q = transitionQ - eBind;
167  G4double cmMomentum = Q*(Q + 2.*daughterMass)/(Q + daughterMass)/2.;
168 
169  G4double costheta = 2.*G4UniformRand() - 1.0;
170  G4double sintheta = std::sqrt(1.0 - costheta*costheta);
171  G4double phi = twopi*G4UniformRand()*rad;
172  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),
173  costheta);
174  G4double KE = cmMomentum;
175  G4DynamicParticle* daughterParticle =
176  new G4DynamicParticle(G4MT_daughters[1], direction, KE, 0.0);
177  products->PushProducts(daughterParticle);
178 
179  KE = std::sqrt(cmMomentum*cmMomentum + daughterMass*daughterMass) - daughterMass;
180  daughterParticle =
181  new G4DynamicParticle(G4MT_daughters[0], -1.0*direction, KE, daughterMass);
182  products->PushProducts(daughterParticle);
183 
184  G4int nArm = armProducts.size();
185  if (nArm > 0) {
186  G4ThreeVector bst = daughterParticle->Get4Momentum().boostVector();
187  for (G4int i = 0; i < nArm; ++i) {
188  G4DynamicParticle* dp = armProducts[i];
189  G4LorentzVector lv = dp->Get4Momentum().boost(bst);
190  dp->Set4Momentum(lv);
191  products->PushProducts(dp);
192  }
193  }
194 
195  // Energy conservation check
196  /*
197  G4int newSize = products->entries();
198  G4DynamicParticle* temp = 0;
199  G4double KEsum = 0.0;
200  for (G4int i = 0; i < newSize; i++) {
201  temp = products->operator[](i);
202  KEsum += temp->GetKineticEnergy();
203  }
204 
205  G4double eCons = (transitionQ - KEsum)/keV;
206  G4cout << " EC check: Ediff (keV) = " << eCons << G4endl;
207  */
208  return products;
209 }
210 
211 
212 void G4ECDecay::DumpNuclearInfo()
213 {
214  G4cout << " G4ECDecay of parent nucleus " << GetParentName() << " from ";
215  if (theMode == 3) {
216  G4cout << "K shell";
217  } else if (theMode == 4) {
218  G4cout << "L shell";
219  } else if (theMode == 5) {
220  G4cout << "M shell";
221  }
222  else if (theMode == 6) {
223  G4cout << "N shell";
224  }
225  G4cout << G4endl;
226  G4cout << " to " << GetDaughterName(0) << " + " << GetDaughterName(1)
227  << " with branching ratio " << GetBR() << "% and Q value "
228  << transitionQ << G4endl;
229 }
230 void G4ECDecay::DefineSubshellProbabilities(G4int Z, G4int )
231 { //Implementation for the case of allowed transitions
232  //PL1+PL2=1. , PM1+PM2=1., PN1+PN2=1.
233  PL1 = 1./(1+PL2overPL1[Z-1]);
234  PL2 = PL1*PL2overPL1[Z-1];
235  PM1 = 1./(1+PM2overPM1[Z-1]);
236  PM2 = PM1*PM2overPM1[Z-1];
237  PN1 = 1./(1+PN2overPN1[Z-1]);
238  PN2 = PN1*PN2overPN1[Z-1];
239 }
241 // Table of subshell ratio probability PL2/PL1 in function of Z
242 // PL2/PL1 = (fL2/gL1)^2
243 // with gL1 and fL2 the bound electron radial wave amplitudes taken from
244 // Bambynek et al., Rev. Modern Physics, vol. 49, 1971, table IX
245 // For Z=18 interpolation between Z=17 and Z=19 to avoid a jump in PL2/Pl1
247 const G4double G4ECDecay::PL2overPL1[100] = {
248 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 1.8722e-04,
249 2.6438e-04, 3.5456e-04, 4.5790e-04, 6.1588e-04, 7.8944e-04, 9.8530e-04, 1.2030e-03,
250 1.4361e-03, 1.6886e-03, 1.9609e-03, 2.2641e-03, 2.5674e-03, 2.9019e-03, 3.2577e-03,
251 3.6338e-03, 4.0310e-03, 4.4541e-03, 4.8943e-03, 5.3620e-03, 5.8523e-03, 6.3650e-03,
252 6.9061e-03, 7.4607e-03, 8.0398e-03, 8.6417e-03, 9.2665e-03, 9.9150e-03, 1.0588e-02,
253 1.1284e-02, 1.2004e-02, 1.2744e-02, 1.3518e-02, 1.4312e-02, 1.5136e-02, 1.5981e-02,
254 1.6857e-02, 1.7764e-02, 1.8696e-02, 1.9682e-02, 2.0642e-02, 2.1661e-02, 2.2708e-02,
255 2.3788e-02, 2.4896e-02, 2.6036e-02, 2.7217e-02, 2.8409e-02, 2.9646e-02, 3.0917e-02,
256 3.2220e-02, 3.3561e-02, 3.4937e-02, 3.6353e-02, 3.7805e-02, 3.9297e-02, 4.0826e-02,
257 4.2399e-02, 4.4010e-02, 4.5668e-02, 4.7368e-02, 4.9115e-02, 5.0896e-02, 5.2744e-02,
258 5.4625e-02, 5.6565e-02, 5.8547e-02, 6.0593e-02, 6.2690e-02, 6.4844e-02, 6.7068e-02,
259 6.9336e-02, 7.1667e-02, 7.4075e-02, 7.6544e-02, 7.9085e-02, 8.1688e-02, 8.4371e-02,
260 8.7135e-02, 8.9995e-02, 9.2919e-02, 9.5949e-02, 9.9036e-02, 1.0226e-01, 1.0555e-01,
261 1.0899e-01, 1.1249e-01, 1.1613e-01, 1.1989e-01, 1.2379e-01, 1.2780e-01, 1.3196e-01,
262 1.3627e-01, 1.4071e-01};
264 // Table of subshell ratio probability PM2/PM1 in function of Z
265 // PM2/PM1 = (fM2/gM1)^2
266 // with gM1 and fM2 the bound electron radial wave amplitudes taken from
267 // Bambynek et al., Rev. Modern Physics, vol. 49, 1971, table IX
269 const G4double G4ECDecay::PM2overPM1[100] = {
270 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
271 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
272 1.0210e-03, 1.2641e-03, 1.5231e-03, 1.7990e-03, 2.1938e-03, 2.5863e-03, 2.9621e-03,
273 3.3637e-03, 3.7909e-03, 4.2049e-03, 4.7021e-03, 5.1791e-03, 5.6766e-03, 6.1952e-03,
274 6.7045e-03, 7.2997e-03, 7.9438e-03, 8.6271e-03, 9.3294e-03, 1.0058e-02, 1.0813e-02,
275 1.1594e-02, 1.2408e-02, 1.3244e-02, 1.4118e-02, 1.5023e-02, 1.5962e-02, 1.6919e-02,
276 1.7910e-02, 1.8934e-02, 1.9986e-02, 2.1072e-02, 2.2186e-02, 2.3336e-02, 2.4524e-02,
277 2.5750e-02, 2.7006e-02, 2.8302e-02, 2.9629e-02, 3.0994e-02, 3.2399e-02, 3.3845e-02,
278 3.5328e-02, 3.6852e-02, 3.8414e-02, 4.0025e-02, 4.1673e-02, 4.3368e-02, 4.5123e-02,
279 4.6909e-02, 4.8767e-02, 5.0662e-02, 5.2612e-02, 5.4612e-02, 5.6662e-02, 5.8773e-02,
280 6.0930e-02, 6.3141e-02, 6.5413e-02, 6.7752e-02, 7.0139e-02, 7.2603e-02, 7.5127e-02,
281 7.7721e-02, 8.0408e-02, 8.3128e-02, 8.5949e-02, 8.8843e-02, 9.1824e-02, 9.4888e-02,
282 9.8025e-02, 1.0130e-01, 1.0463e-01, 1.0806e-01, 1.1159e-01, 1.1526e-01, 1.1900e-01,
283 1.2290e-01, 1.2688e-01, 1.3101e-01, 1.3528e-01, 1.3969e-01, 1.4425e-01, 1.4896e-01,
284 1.5384e-01, 1.5887e-01};
286 // Table of subshell ratio probability PN2/PN1 in function of Z
287 // PN2/PN1 = (fN2/gN1)^2
288 // with gN1 and fN2 are the bound electron radial wave amplitude taken from
289 // Bambynek et al., Rev. Modern Physics, vol. 49, 1971, table IX
290 // For Z=44 interpolation between Z=43 and Z=45 to avoid a jump in PN2/PN1
292 const G4double G4ECDecay::PN2overPN1[100] = {
293 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
294 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
295 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
296 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
297 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00, 0.0000e+00,
298 0.0000e+00, 9.6988e-03, 1.0797e-02, 1.1706e-02, 1.2603e-02, 1.3408e-02, 1.4352e-02,
299 1.5511e-02, 1.6579e-02, 1.7646e-02, 1.8731e-02, 1.9886e-02, 2.1069e-02, 2.2359e-02,
300 2.3710e-02, 2.5058e-02, 2.6438e-02, 2.7843e-02, 2.9283e-02, 3.0762e-02, 3.2275e-02,
301 3.3843e-02, 3.5377e-02, 3.6886e-02, 3.8502e-02, 4.0159e-02, 4.1867e-02, 4.3617e-02,
302 4.5470e-02, 4.7247e-02, 4.9138e-02, 5.1065e-02, 5.3049e-02, 5.5085e-02, 5.7173e-02,
303 5.9366e-02, 6.1800e-02, 6.3945e-02, 6.6333e-02, 6.8785e-02, 7.1303e-02, 7.3801e-02,
304 7.6538e-02, 7.9276e-02, 8.2070e-02, 8.4959e-02, 8.7940e-02, 9.0990e-02, 9.4124e-02,
305 9.7337e-02, 1.0069e-01, 1.0410e-01, 1.0761e-01, 1.1122e-01, 1.1499e-01, 1.1882e-01,
306 1.2282e-01, 1.2709e-01, 1.3114e-01, 1.3549e-01, 1.4001e-01, 1.4465e-01, 1.4946e-01,
307 1.5443e-01, 1.5954e-01};
308