| 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 // File: G4BetaMinusDecay.cc
29 // Author: D.H. Wright (SLAC)
30 // Date: 25 October 2014
31 // Modifications:
32 // 23.08.2023 V.Ivanchenko make it thread s
33 //
34 //////////////////////////////////////////////
35
36 #include "G4BetaMinusDecay.hh"
37 #include "G4BetaDecayCorrections.hh"
38 #include "G4ThreeVector.hh"
39 #include "G4LorentzVector.hh"
40 #include "G4DynamicParticle.hh"
41 #include "G4DecayProducts.hh"
42 #include "G4PhysicalConstants.hh"
43 #include "G4SystemOfUnits.hh"
44 #include "G4Electron.hh"
45 #include "G4AntiNeutrinoE.hh"
46 #include "G4RandomDirection.hh"
47 #include "G4BetaSpectrumSampler.hh"
48 #include <iostream>
49 #include <iomanip>
50
51 namespace {
52 const G4double eMass = CLHEP::electron_mass_
53 }
54
55 G4BetaMinusDecay::G4BetaMinusDecay(const G4Par
56 const G4dou
57 const G4dou
58 const G4Ion
59 const G4Bet
60 : G4NuclearDecay("beta- decay", BetaMinus, ex
61 maxEnergy(e0/eMass),
62 estep(maxEnergy/(G4double)(npti - 1))
63 {
64 SetParent(theParentNucleus); // Store name
65 SetBR(branch);
66 SetNumberOfDaughters(3);
67
68 fPrimaryIon = theParentNucleus;
69 fLepton = G4Electron::Electron();
70 fNeutrino = G4AntiNeutrinoE::AntiNeutrinoE()
71
72 G4IonTable* theIonTable = G4ParticleTable::G
73 G4int daughterZ = theParentNucleus->GetAtomi
74 G4int daughterA = theParentNucleus->GetAtomi
75 fResIon = const_cast<const G4ParticleDefinit
76
77 parentMass = theParentNucleus->GetPDGMass();
78 resMass = fResIon->GetPDGMass();
79
80 SetUpBetaSpectrumSampler(daughterZ, daughter
81
82 SetDaughter(0, fResIon);
83 SetDaughter(1, fLepton);
84 SetDaughter(2, fNeutrino);
85
86 // Fill G4MT_parent with theParentNucleus (s
87 CheckAndFillParent();
88
89 // Fill G4MT_daughters with e-, nu and resid
90 CheckAndFillDaughters();
91 }
92
93 G4DecayProducts* G4BetaMinusDecay::DecayIt(G4d
94 {
95 // Set up final state
96 // parentParticle is set at rest here becaus
97 // is done later
98 G4DynamicParticle prim(fPrimaryIon, G4ThreeV
99 G4DecayProducts* products = new G4DecayProdu
100
101 // Generate electron isotropic in angle, wit
102 const G4double eKE = eMass*G4BetaSpectrumSam
103
104 G4double eMomentum = std::sqrt(eKE*(eKE + 2.
105 G4ThreeVector dir = G4RandomDirection();
106 G4DynamicParticle* dp = new G4DynamicParticl
107 products->PushProducts(dp);
108 /*
109 G4cout << "G4BetaPlusDecay::DecayIt: " << fP
110 << " -> " << fResIon->GetParticleName() <<
111 << " + " << fNeutrino->GetParticleName() <<
112 << G4endl;
113 */
114
115 // 4-momentum of residual ion and neutrino
116 G4LorentzVector lv(-eMomentum*dir.x(), -eMom
117 parentMass - eKE - eMass)
118
119 // centrum of mass system
120 G4double M = lv.mag();
121 const G4double elim = CLHEP::eV;
122 G4double edel = M - resMass;
123 // Free energy should be above limit
124 if (edel >= elim) {
125 // neutrino
126 G4double eNu = 0.5*(M - resMass*resMass/M)
127 G4LorentzVector lvnu(eNu*G4RandomDirection
128 lvnu.boost(lv.boostVector());
129 dir = lvnu.vect().unit();
130 dp = new G4DynamicParticle(fNeutrino, dir,
131 products->PushProducts(dp);
132
133 // residual
134 lv -= lvnu;
135 dir = lv.vect().unit();
136 G4double ekin = std::max(lv.e() - resMass,
137 dp = new G4DynamicParticle(fResIon, dir, e
138 products->PushProducts(dp);
139
140 } else {
141 // neglecting relativistic kinematic and g
142 dp = new G4DynamicParticle(fNeutrino, G4Ra
143 products->PushProducts(dp);
144 dp = new G4DynamicParticle(fResIon, G4Thre
145 products->PushProducts(dp);
146 }
147
148 return products;
149 }
150
151
152 void
153 G4BetaMinusDecay::SetUpBetaSpectrumSampler(con
154 con
155 con
156 {
157 cdf[0] = 0.0;
158
159 // Check for cases in which Q < 0
160 if (maxEnergy > 0.) {
161 G4BetaDecayCorrections corrections(daughte
162
163 // Fill array to store cumulative spectrum
164 G4double ex; // Kinetic energy normalized
165 G4double p; // Momentum in units of elec
166 G4double f; // Spectral shape function
167 G4double f0 = 0.0;
168 G4double sum = 0.0;
169 for (G4int i = 1; i < npti-1; ++i) {
170 ex = estep*i;
171 p = std::sqrt(ex*(ex + 2.));
172 f = p*(1. + ex)*(maxEnergy - ex)*(maxEne
173
174 // Apply Fermi factor to get allowed sha
175 f *= corrections.FermiFunction(1. + ex);
176
177 // Apply shape factor for forbidden tran
178 f *= corrections.ShapeFactor(betaType, p
179 sum += f + f0;
180 cdf[i] = sum;
181 f0 = f;
182 }
183 cdf[npti-1] = sum + f0;
184 } else {
185 for (G4int i = 1; i < npti; ++i) { cdf[i]
186 }
187 }
188
189
190 void G4BetaMinusDecay::DumpNuclearInfo()
191 {
192 G4cout << " G4BetaMinusDecay " << fPrimaryI
193 << " -> " << fResIon->GetParticleName() <<
194 << " + " << fNeutrino->GetParticleName() <<
195 << maxEnergy*eMass << " BR=" << GetBR() <<
196 }
197
198