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: G4LFission.cc,v 1.15 2007/02/26 19:29:30 dennis Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-02-ref-02 $ |
---|
29 | // |
---|
30 | // |
---|
31 | // G4 Model: Low Energy Fission |
---|
32 | // F.W. Jones, TRIUMF, 03-DEC-96 |
---|
33 | // |
---|
34 | // This is a prototype of a low-energy fission process. |
---|
35 | // Currently it is based on the GHEISHA routine FISSIO, |
---|
36 | // and conforms fairly closely to the original Fortran. |
---|
37 | // Note: energy is in MeV and momentum is in MeV/c. |
---|
38 | // |
---|
39 | // use -scheme for elastic scattering: HPW, 20th June 1997 |
---|
40 | // the code comes mostly from the old Low-energy Fission class |
---|
41 | // |
---|
42 | // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange. |
---|
43 | // |
---|
44 | |
---|
45 | #include "globals.hh" |
---|
46 | #include "G4LFission.hh" |
---|
47 | #include "Randomize.hh" |
---|
48 | |
---|
49 | G4LFission::G4LFission() : G4HadronicInteraction("G4LFission") |
---|
50 | { |
---|
51 | init(); |
---|
52 | SetMinEnergy( 0.0*GeV ); |
---|
53 | SetMaxEnergy( DBL_MAX ); |
---|
54 | |
---|
55 | } |
---|
56 | |
---|
57 | G4LFission::~G4LFission() |
---|
58 | { |
---|
59 | theParticleChange.Clear(); |
---|
60 | } |
---|
61 | |
---|
62 | |
---|
63 | void |
---|
64 | G4LFission::init() |
---|
65 | { |
---|
66 | G4int i; |
---|
67 | G4double xx = 1. - 0.5; |
---|
68 | G4double xxx = std::sqrt(2.29*xx); |
---|
69 | spneut[0] = std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.; |
---|
70 | for (i = 2; i <= 10; i++) { |
---|
71 | xx = i*1. - 0.5; |
---|
72 | xxx = std::sqrt(2.29*xx); |
---|
73 | spneut[i-1] = spneut[i-2] + std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.; |
---|
74 | } |
---|
75 | for (i = 1; i <= 10; i++) { |
---|
76 | spneut[i-1] = spneut[i-1]/spneut[9]; |
---|
77 | if (verboseLevel > 1) G4cout << "G4LFission::init: i=" << i << |
---|
78 | " spneut=" << spneut[i-1] << G4endl; |
---|
79 | } |
---|
80 | } |
---|
81 | |
---|
82 | |
---|
83 | G4HadFinalState* |
---|
84 | G4LFission::ApplyYourself(const G4HadProjectile & aTrack,G4Nucleus & targetNucleus) |
---|
85 | { |
---|
86 | theParticleChange.Clear(); |
---|
87 | const G4HadProjectile* aParticle = &aTrack; |
---|
88 | |
---|
89 | G4double N = targetNucleus.GetN(); |
---|
90 | G4double Z = targetNucleus.GetZ(); |
---|
91 | theParticleChange.SetStatusChange(stopAndKill); |
---|
92 | |
---|
93 | G4double P = aParticle->GetTotalMomentum()/MeV; |
---|
94 | G4double Px = aParticle->Get4Momentum().vect().x(); |
---|
95 | G4double Py = aParticle->Get4Momentum().vect().y(); |
---|
96 | G4double Pz = aParticle->Get4Momentum().vect().z(); |
---|
97 | G4double E = aParticle->GetTotalEnergy()/MeV; |
---|
98 | G4double E0 = aParticle->GetDefinition()->GetPDGMass()/MeV; |
---|
99 | G4double Q = aParticle->GetDefinition()->GetPDGCharge(); |
---|
100 | if (verboseLevel > 1) { |
---|
101 | G4cout << "G4LFission:ApplyYourself: incident particle:" << G4endl; |
---|
102 | G4cout << "P " << P << " MeV/c" << G4endl; |
---|
103 | G4cout << "Px " << Px << " MeV/c" << G4endl; |
---|
104 | G4cout << "Py " << Py << " MeV/c" << G4endl; |
---|
105 | G4cout << "Pz " << Pz << " MeV/c" << G4endl; |
---|
106 | G4cout << "E " << E << " MeV" << G4endl; |
---|
107 | G4cout << "mass " << E0 << " MeV" << G4endl; |
---|
108 | G4cout << "charge " << Q << G4endl; |
---|
109 | } |
---|
110 | // GHEISHA ADD operation to get total energy, mass, charge: |
---|
111 | if (verboseLevel > 1) { |
---|
112 | G4cout << "G4LFission:ApplyYourself: material:" << G4endl; |
---|
113 | G4cout << "A " << N << G4endl; |
---|
114 | G4cout << "Z " << Z << G4endl; |
---|
115 | G4cout << "atomic mass " << |
---|
116 | Atomas(N, Z) << "MeV" << G4endl; |
---|
117 | } |
---|
118 | E = E + Atomas(N, Z); |
---|
119 | G4double E02 = E*E - P*P; |
---|
120 | E0 = std::sqrt(std::abs(E02)); |
---|
121 | if (E02 < 0) E0 = -E0; |
---|
122 | Q = Q + Z; |
---|
123 | if (verboseLevel > 1) { |
---|
124 | G4cout << "G4LFission:ApplyYourself: total:" << G4endl; |
---|
125 | G4cout << "E " << E << " MeV" << G4endl; |
---|
126 | G4cout << "mass " << E0 << " MeV" << G4endl; |
---|
127 | G4cout << "charge " << Q << G4endl; |
---|
128 | } |
---|
129 | Px = -Px; |
---|
130 | Py = -Py; |
---|
131 | Pz = -Pz; |
---|
132 | |
---|
133 | G4double e1 = aParticle->GetKineticEnergy()/MeV; |
---|
134 | if (e1 < 1.) e1 = 1.; |
---|
135 | |
---|
136 | // Average number of neutrons |
---|
137 | G4double avern = 2.569 + 0.559*std::log(e1); |
---|
138 | G4bool photofission = 0; // For now |
---|
139 | // Take the following value if photofission is not included |
---|
140 | if (!photofission) avern = 2.569 + 0.900*std::log(e1); |
---|
141 | |
---|
142 | // Average number of gammas |
---|
143 | G4double averg = 9.500 + 0.600*std::log(e1); |
---|
144 | |
---|
145 | G4double ran = G4RandGauss::shoot(); |
---|
146 | // Number of neutrons |
---|
147 | G4int nn = static_cast<G4int>(avern + ran*1.23 + 0.5); |
---|
148 | ran = G4RandGauss::shoot(); |
---|
149 | // Number of gammas |
---|
150 | G4int ng = static_cast<G4int>(averg + ran*3. + 0.5); |
---|
151 | if (nn < 1) nn = 1; |
---|
152 | if (ng < 1) ng = 1; |
---|
153 | G4double exn = 0.; |
---|
154 | G4double exg = 0.; |
---|
155 | |
---|
156 | // Make secondary neutrons and distribute kinetic energy |
---|
157 | G4DynamicParticle* aNeutron; |
---|
158 | G4int i; |
---|
159 | for (i = 1; i <= nn; i++) { |
---|
160 | ran = G4UniformRand(); |
---|
161 | G4int j; |
---|
162 | for (j = 1; j <= 10; j++) { |
---|
163 | if (ran < spneut[j-1]) goto label12; |
---|
164 | } |
---|
165 | j = 10; |
---|
166 | label12: |
---|
167 | ran = G4UniformRand(); |
---|
168 | G4double ekin = (j - 1)*1. + ran; |
---|
169 | exn = exn + ekin; |
---|
170 | aNeutron = new G4DynamicParticle(G4Neutron::NeutronDefinition(), |
---|
171 | G4ParticleMomentum(1.,0.,0.), |
---|
172 | ekin*MeV); |
---|
173 | theParticleChange.AddSecondary(aNeutron); |
---|
174 | } |
---|
175 | |
---|
176 | // Make secondary gammas and distribute kinetic energy |
---|
177 | G4DynamicParticle* aGamma; |
---|
178 | for (i = 1; i <= ng; i++) { |
---|
179 | ran = G4UniformRand(); |
---|
180 | G4double ekin = -0.87*std::log(ran); |
---|
181 | exg = exg + ekin; |
---|
182 | aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(), |
---|
183 | G4ParticleMomentum(1.,0.,0.), |
---|
184 | ekin*MeV); |
---|
185 | theParticleChange.AddSecondary(aGamma); |
---|
186 | } |
---|
187 | |
---|
188 | // Distribute momentum vectors and do Lorentz transformation |
---|
189 | |
---|
190 | G4HadSecondary* theSecondary; |
---|
191 | |
---|
192 | for (i = 1; i <= nn + ng; i++) { |
---|
193 | G4double ran1 = G4UniformRand(); |
---|
194 | G4double ran2 = G4UniformRand(); |
---|
195 | G4double cost = -1. + 2.*ran1; |
---|
196 | G4double sint = std::sqrt(std::abs(1. - cost*cost)); |
---|
197 | G4double phi = ran2*twopi; |
---|
198 | // G4cout << ran1 << " " << ran2 << G4endl; |
---|
199 | // G4cout << cost << " " << sint << " " << phi << G4endl; |
---|
200 | theSecondary = theParticleChange.GetSecondary(i - 1); |
---|
201 | G4double pp = theSecondary->GetParticle()->GetTotalMomentum()/MeV; |
---|
202 | G4double px = pp*sint*std::sin(phi); |
---|
203 | G4double py = pp*sint*std::cos(phi); |
---|
204 | G4double pz = pp*cost; |
---|
205 | // G4cout << pp << G4endl; |
---|
206 | // G4cout << px << " " << py << " " << pz << G4endl; |
---|
207 | G4double e = theSecondary->GetParticle()->GetTotalEnergy()/MeV; |
---|
208 | G4double e0 = theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV; |
---|
209 | |
---|
210 | G4double a = px*Px + py*Py + pz*Pz; |
---|
211 | a = (a/(E + E0) - e)/E0; |
---|
212 | |
---|
213 | px = px + a*Px; |
---|
214 | py = py + a*Py; |
---|
215 | pz = pz + a*Pz; |
---|
216 | G4double p2 = px*px + py*py + pz*pz; |
---|
217 | pp = std::sqrt(p2); |
---|
218 | e = std::sqrt(e0*e0 + p2); |
---|
219 | G4double ekin = e - theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV; |
---|
220 | theSecondary->GetParticle()->SetMomentumDirection(G4ParticleMomentum(px/pp, |
---|
221 | py/pp, |
---|
222 | pz/pp)); |
---|
223 | theSecondary->GetParticle()->SetKineticEnergy(ekin*MeV); |
---|
224 | } |
---|
225 | |
---|
226 | return &theParticleChange; |
---|
227 | } |
---|
228 | |
---|
229 | // Computes atomic mass in MeV (translation of GHEISHA routine ATOMAS) |
---|
230 | // Not optimized: conforms closely to original Fortran. |
---|
231 | |
---|
232 | G4double |
---|
233 | G4LFission::Atomas(const G4double A, const G4double Z) |
---|
234 | { |
---|
235 | G4double rmel = G4Electron::ElectronDefinition()->GetPDGMass()/MeV; |
---|
236 | G4double rmp = G4Proton::ProtonDefinition()->GetPDGMass()/MeV; |
---|
237 | G4double rmn = G4Neutron::NeutronDefinition()->GetPDGMass()/MeV; |
---|
238 | G4double rmd = G4Deuteron::DeuteronDefinition()->GetPDGMass()/MeV; |
---|
239 | G4double rma = G4Alpha::AlphaDefinition()->GetPDGMass()/MeV; |
---|
240 | |
---|
241 | G4int ia = static_cast<G4int>(A + 0.5); |
---|
242 | if (ia < 1) return 0; |
---|
243 | G4int iz = static_cast<G4int>(Z + 0.5); |
---|
244 | if (iz < 0) return 0; |
---|
245 | if (iz > ia) return 0; |
---|
246 | |
---|
247 | if (ia == 1) { |
---|
248 | if (iz == 0) return rmn; //neutron |
---|
249 | if (iz == 1) return rmp + rmel; //Hydrogen |
---|
250 | } |
---|
251 | else if (ia == 2 && iz == 1) { |
---|
252 | return rmd; //Deuteron |
---|
253 | } |
---|
254 | else if (ia == 4 && iz == 2) { |
---|
255 | return rma; //Alpha |
---|
256 | } |
---|
257 | |
---|
258 | G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel |
---|
259 | - 15.67*A |
---|
260 | + 17.23*std::pow(A, 2./3.) |
---|
261 | + 93.15*(A/2. - Z)*(A/2. - Z)/A |
---|
262 | + 0.6984523*Z*Z/std::pow(A, 1./3.); |
---|
263 | G4int ipp = (ia - iz)%2; |
---|
264 | G4int izz = iz%2; |
---|
265 | if (ipp == izz) mass = mass + (ipp + izz -1)*12.*std::pow(A, -0.5); |
---|
266 | |
---|
267 | return mass; |
---|
268 | } |
---|