Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ECDecay.cc
Go to the documentation of this file.
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 
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 
59  G4IonTable* theIonTable =
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 }
66 
67 
69 {}
70 
71 
73 {
74  // Fill G4MT_parent with theParentNucleus (stored by SetParent in ctor)
76 
77  // Fill G4MT_daughters with alpha and residual nucleus (stored by SetDaughter)
79 
80  // Get shell number of captured electron
81  G4int shellIndex = -1;
82  switch (theMode)
83  {
84  case KshellEC:
85  shellIndex = 0;
86  break;
87  case LshellEC:
88  shellIndex = G4int(G4UniformRand()*3) + 1;
89  break;
90  case MshellEC:
91  shellIndex = G4int(G4UniformRand()*3) + 4;
92  break;
93  default:
94  G4Exception("G4ECDecay::DecayIt()", "HAD_RDM_009",
95  FatalException, "Invalid electron shell selected");
96  }
97 
98  // Initialize decay products with parent nucleus at rest
99  G4DynamicParticle parentParticle(G4MT_parent, G4ThreeVector(0,0,0), 0.0);
100  G4DecayProducts* products = new G4DecayProducts(parentParticle);
101  G4double eBind = 0.0;
102 
103  // G4LossTableManager must already be initialized with G4UAtomicDeexcitation
104  // This is currently done in G4RadioactiveDecay::BuildPhysicsTable
105  G4VAtomDeexcitation* atomDeex =
107  std::vector<G4DynamicParticle*> armProducts;
108 
109  if (applyARM) {
110  if (atomDeex) {
113  if (shellIndex >= nShells) shellIndex = nShells;
115  const G4AtomicShell* shell = atomDeex->GetAtomicShell(aZ, as);
116  eBind = shell->BindingEnergy();
117  if (atomDeex->IsFluoActive() && aZ > 5 && aZ < 100) {
118  // Do atomic relaxation
119  // VI, SI
120  // Allows fixing of Bugzilla 1727
121  //const G4double deexLimit = 0.1*keV;
122  G4double deexLimit = 0.1*keV;
123  if (G4EmParameters::Instance()->DeexcitationIgnoreCut()) deexLimit =0.;
124  //
125  atomDeex->GenerateParticles(&armProducts, shell, aZ, deexLimit, deexLimit);
126  }
127 
128  G4double productEnergy = 0.;
129  for (G4int i = 0; i < G4int(armProducts.size()); i++)
130  productEnergy += armProducts[i]->GetKineticEnergy();
131 
132  G4double deficit = shell->BindingEnergy() - productEnergy;
133  if (deficit > 0.0) {
134  // Add a dummy electron to make up extra energy
135  G4double cosTh = 1.-2.*G4UniformRand();
136  G4double sinTh = std::sqrt(1.- cosTh*cosTh);
137  G4double phi = twopi*G4UniformRand();
138 
139  G4ThreeVector electronDirection(sinTh*std::sin(phi),
140  sinTh*std::cos(phi), cosTh);
141  G4DynamicParticle* extra =
142  new G4DynamicParticle(G4Electron::Electron(), electronDirection,
143  deficit);
144  armProducts.push_back(extra);
145  }
146  } // atomDeex
147  } // applyARM
148 
149  G4double daughterMass = G4MT_daughters[0]->GetPDGMass();
150 
151  // CM momentum using Q value corrected for binding energy of captured electron
152  G4double Q = transitionQ - eBind;
153  G4double cmMomentum = Q*(Q + 2.*daughterMass)/(Q + daughterMass)/2.;
154 
155  G4double costheta = 2.*G4UniformRand() - 1.0;
156  G4double sintheta = std::sqrt(1.0 - costheta*costheta);
158  G4ThreeVector direction(sintheta*std::cos(phi),sintheta*std::sin(phi),
159  costheta);
160  G4double KE = cmMomentum;
161  G4DynamicParticle* daughterParticle =
162  new G4DynamicParticle(G4MT_daughters[1], direction, KE, 0.0);
163  products->PushProducts(daughterParticle);
164 
165  KE = std::sqrt(cmMomentum*cmMomentum + daughterMass*daughterMass) - daughterMass;
166  daughterParticle =
167  new G4DynamicParticle(G4MT_daughters[0], -1.0*direction, KE, daughterMass);
168  products->PushProducts(daughterParticle);
169 
170  G4int nArm = armProducts.size();
171  if (nArm > 0) {
172  G4ThreeVector bst = daughterParticle->Get4Momentum().boostVector();
173  for (G4int i = 0; i < nArm; ++i) {
174  G4DynamicParticle* dp = armProducts[i];
175  G4LorentzVector lv = dp->Get4Momentum().boost(bst);
176  dp->Set4Momentum(lv);
177  products->PushProducts(dp);
178  }
179  }
180 
181  // Energy conservation check
182  /*
183  G4int newSize = products->entries();
184  G4DynamicParticle* temp = 0;
185  G4double KEsum = 0.0;
186  for (G4int i = 0; i < newSize; i++) {
187  temp = products->operator[](i);
188  KEsum += temp->GetKineticEnergy();
189  }
190 
191  G4double eCons = (transitionQ - KEsum)/keV;
192  G4cout << " EC check: Ediff (keV) = " << eCons << G4endl;
193  */
194  return products;
195 }
196 
197 
199 {
200  G4cout << " G4ECDecay of parent nucleus " << GetParentName() << " from ";
201  if (theMode == 3) {
202  G4cout << "K shell";
203  } else if (theMode == 4) {
204  G4cout << "L shell";
205  } else if (theMode == 5) {
206  G4cout << "M shell";
207  }
208  G4cout << G4endl;
209  G4cout << " to " << GetDaughterName(0) << " + " << GetDaughterName(1)
210  << " with branching ratio " << GetBR() << "% and Q value "
211  << transitionQ << G4endl;
212 }
213 
const G4RadioactiveDecayMode theMode
void CheckAndFillDaughters()
G4double GetBR() const
Hep3Vector boostVector() const
void SetBR(G4double value)
static G4LossTableManager * Instance()
G4bool IsFluoActive() const
CLHEP::Hep3Vector G4ThreeVector
const G4String & GetParentName() const
G4int PushProducts(G4DynamicParticle *aParticle)
G4ParticleDefinition * G4MT_parent
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:503
G4ParticleDefinition ** G4MT_daughters
static double Q[]
static constexpr double rad
Definition: G4SIunits.hh:149
G4double BindingEnergy() const
int G4int
Definition: G4Types.hh:78
G4bool DeexcitationIgnoreCut() const
G4int GetAtomicNumber() const
static constexpr double twopi
Definition: G4SIunits.hh:76
virtual const G4AtomicShell * GetAtomicShell(G4int Z, G4AtomicShellEnumerator shell)=0
G4IonTable * GetIonTable() const
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
void SetNumberOfDaughters(G4int value)
HepLorentzVector & boost(double, double, double)
G4RadioactiveDecayMode
G4int GetAtomicMass() const
G4LorentzVector Get4Momentum() const
const G4String & GetDaughterName(G4int anIndex) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
void Set4Momentum(const G4LorentzVector &momentum)
G4ECDecay(const G4ParticleDefinition *theParentNucleus, const G4double &theBR, const G4double &Qvalue, const G4double &excitation, const G4Ions::G4FloatLevelBase &flb, const G4RadioactiveDecayMode &mode)
Definition: G4ECDecay.cc:47
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
void SetParent(const G4ParticleDefinition *particle_type)
static G4EmParameters * Instance()
void SetDaughter(G4int anIndex, const G4ParticleDefinition *particle_type)
virtual void DumpNuclearInfo()
Definition: G4ECDecay.cc:198
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
G4VAtomDeexcitation * AtomDeexcitation()
double G4double
Definition: G4Types.hh:76
void GenerateParticles(std::vector< G4DynamicParticle * > *secVect, const G4AtomicShell *, G4int Z, G4int coupleIndex)
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4ECDecay.cc:72
static constexpr double keV
Definition: G4SIunits.hh:216
virtual ~G4ECDecay()
Definition: G4ECDecay.cc:68
G4AtomicShellEnumerator
G4FloatLevelBase
Definition: G4Ions.hh:95
static G4int GetNumberOfShells(G4int Z)