| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // G4 Model: Low Energy Fission 28 // G4 Model: Low Energy Fission
29 // F.W. Jones, TRIUMF, 03-DEC-96 29 // F.W. Jones, TRIUMF, 03-DEC-96
30 // 30 //
31 // This is a prototype of a low-energy fission 31 // This is a prototype of a low-energy fission process.
32 // Currently it is based on the GHEISHA routin 32 // Currently it is based on the GHEISHA routine FISSIO,
33 // and conforms fairly closely to the original 33 // and conforms fairly closely to the original Fortran.
34 // Note: energy is in MeV and momentum is in M 34 // Note: energy is in MeV and momentum is in MeV/c.
35 // 35 //
36 // use -scheme for elastic scattering: HPW, 20 36 // use -scheme for elastic scattering: HPW, 20th June 1997
37 // the code comes mostly from the old Low-ener 37 // the code comes mostly from the old Low-energy Fission class
38 // 38 //
39 // 25-JUN-98 FWJ: replaced missing Initialize 39 // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
40 40
41 #include <iostream> 41 #include <iostream>
42 42
43 #include "G4LFission.hh" 43 #include "G4LFission.hh"
44 #include "globals.hh" 44 #include "globals.hh"
45 #include "G4Exp.hh" 45 #include "G4Exp.hh"
46 #include "G4Log.hh" 46 #include "G4Log.hh"
47 #include "G4Pow.hh" 47 #include "G4Pow.hh"
48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
50 #include "Randomize.hh" 50 #include "Randomize.hh"
51 #include "G4PhysicsModelCatalog.hh" <<
52 51
53 G4LFission::G4LFission(const G4String& name) 52 G4LFission::G4LFission(const G4String& name)
54 : G4HadronicInteraction(name), secID(-1) << 53 : G4HadronicInteraction(name)
55 { 54 {
56 init(); 55 init();
57 SetMinEnergy(0.0*GeV); 56 SetMinEnergy(0.0*GeV);
58 SetMaxEnergy(DBL_MAX); 57 SetMaxEnergy(DBL_MAX);
59 G4PhysicsModelCatalog::GetModelID( "model_" <<
60 } 58 }
61 59
62 60
63 G4LFission::~G4LFission() 61 G4LFission::~G4LFission()
64 { 62 {
65 theParticleChange.Clear(); 63 theParticleChange.Clear();
66 } 64 }
67 65
68 66
69 void G4LFission::ModelDescription(std::ostream 67 void G4LFission::ModelDescription(std::ostream& outFile) const
70 { 68 {
71 outFile << "G4LFission is one of the Low Ene 69 outFile << "G4LFission is one of the Low Energy Parameterized\n"
72 << "(LEP) models used to implement n 70 << "(LEP) models used to implement neutron-induced fission of\n"
73 << "nuclei. It is a re-engineered v 71 << "nuclei. It is a re-engineered version of the GHEISHA code\n"
74 << "of H. Fesefeldt which emits neut 72 << "of H. Fesefeldt which emits neutrons and gammas but no\n"
75 << "nuclear fragments. The model is 73 << "nuclear fragments. The model is applicable to all incident\n"
76 << "neutron energies.\n"; 74 << "neutron energies.\n";
77 } 75 }
78 76
79 void G4LFission::init() 77 void G4LFission::init()
80 { 78 {
81 G4int i; 79 G4int i;
82 G4double xx = 1. - 0.5; 80 G4double xx = 1. - 0.5;
83 G4double xxx = std::sqrt(2.29*xx); 81 G4double xxx = std::sqrt(2.29*xx);
84 spneut[0] = G4Exp(-xx/0.965)*(G4Exp(xxx) - 82 spneut[0] = G4Exp(-xx/0.965)*(G4Exp(xxx) - G4Exp(-xxx))/2.;
85 for (i = 2; i <= 10; i++) { 83 for (i = 2; i <= 10; i++) {
86 xx = i*1. - 0.5; 84 xx = i*1. - 0.5;
87 xxx = std::sqrt(2.29*xx); 85 xxx = std::sqrt(2.29*xx);
88 spneut[i-1] = spneut[i-2] + G4Exp(-xx/0. 86 spneut[i-1] = spneut[i-2] + G4Exp(-xx/0.965)*(G4Exp(xxx) - G4Exp(-xxx))/2.;
89 } 87 }
90 for (i = 1; i <= 10; i++) { 88 for (i = 1; i <= 10; i++) {
91 spneut[i-1] = spneut[i-1]/spneut[9]; 89 spneut[i-1] = spneut[i-1]/spneut[9];
92 if (verboseLevel > 1) G4cout << "G4LFiss 90 if (verboseLevel > 1) G4cout << "G4LFission::init: i=" << i <<
93 " spneut=" << spneut[i-1] << G4endl; 91 " spneut=" << spneut[i-1] << G4endl;
94 } 92 }
95 } 93 }
96 94
97 95
98 G4HadFinalState* G4LFission::ApplyYourself(con 96 G4HadFinalState* G4LFission::ApplyYourself(const G4HadProjectile& aTrack,
99 G4N 97 G4Nucleus& targetNucleus)
100 { 98 {
101 theParticleChange.Clear(); 99 theParticleChange.Clear();
102 const G4HadProjectile* aParticle = &aTrack; 100 const G4HadProjectile* aParticle = &aTrack;
103 101
104 G4double N = targetNucleus.GetA_asInt(); 102 G4double N = targetNucleus.GetA_asInt();
105 G4double Z = targetNucleus.GetZ_asInt(); 103 G4double Z = targetNucleus.GetZ_asInt();
106 theParticleChange.SetStatusChange(stopAndKil 104 theParticleChange.SetStatusChange(stopAndKill);
107 105
108 G4double P = aParticle->GetTotalMomentum()/M 106 G4double P = aParticle->GetTotalMomentum()/MeV;
109 G4double Px = aParticle->Get4Momentum().vect 107 G4double Px = aParticle->Get4Momentum().vect().x();
110 G4double Py = aParticle->Get4Momentum().vect 108 G4double Py = aParticle->Get4Momentum().vect().y();
111 G4double Pz = aParticle->Get4Momentum().vect 109 G4double Pz = aParticle->Get4Momentum().vect().z();
112 G4double E = aParticle->GetTotalEnergy()/MeV 110 G4double E = aParticle->GetTotalEnergy()/MeV;
113 G4double E0 = aParticle->GetDefinition()->Ge 111 G4double E0 = aParticle->GetDefinition()->GetPDGMass()/MeV;
114 G4double Q = aParticle->GetDefinition()->Get 112 G4double Q = aParticle->GetDefinition()->GetPDGCharge();
115 if (verboseLevel > 1) { 113 if (verboseLevel > 1) {
116 G4cout << "G4LFission:ApplyYourself: inc 114 G4cout << "G4LFission:ApplyYourself: incident particle:" << G4endl;
117 G4cout << "P " << P << " MeV/c" << 115 G4cout << "P " << P << " MeV/c" << G4endl;
118 G4cout << "Px " << Px << " MeV/c" << 116 G4cout << "Px " << Px << " MeV/c" << G4endl;
119 G4cout << "Py " << Py << " MeV/c" << 117 G4cout << "Py " << Py << " MeV/c" << G4endl;
120 G4cout << "Pz " << Pz << " MeV/c" << 118 G4cout << "Pz " << Pz << " MeV/c" << G4endl;
121 G4cout << "E " << E << " MeV" << G4 119 G4cout << "E " << E << " MeV" << G4endl;
122 G4cout << "mass " << E0 << " MeV" << G 120 G4cout << "mass " << E0 << " MeV" << G4endl;
123 G4cout << "charge " << Q << G4endl; 121 G4cout << "charge " << Q << G4endl;
124 } 122 }
125 // GHEISHA ADD operation to get total energy 123 // GHEISHA ADD operation to get total energy, mass, charge:
126 if (verboseLevel > 1) { 124 if (verboseLevel > 1) {
127 G4cout << "G4LFission:ApplyYourself: mat 125 G4cout << "G4LFission:ApplyYourself: material:" << G4endl;
128 G4cout << "A " << N << G4endl; 126 G4cout << "A " << N << G4endl;
129 G4cout << "Z " << Z << G4endl; 127 G4cout << "Z " << Z << G4endl;
130 G4cout << "atomic mass " << 128 G4cout << "atomic mass " <<
131 Atomas(N, Z) << "MeV" << G4endl; 129 Atomas(N, Z) << "MeV" << G4endl;
132 } 130 }
133 E = E + Atomas(N, Z); 131 E = E + Atomas(N, Z);
134 G4double E02 = E*E - P*P; 132 G4double E02 = E*E - P*P;
135 E0 = std::sqrt(std::abs(E02)); 133 E0 = std::sqrt(std::abs(E02));
136 if (E02 < 0) E0 = -E0; 134 if (E02 < 0) E0 = -E0;
137 Q = Q + Z; 135 Q = Q + Z;
138 if (verboseLevel > 1) { 136 if (verboseLevel > 1) {
139 G4cout << "G4LFission:ApplyYourself: tot 137 G4cout << "G4LFission:ApplyYourself: total:" << G4endl;
140 G4cout << "E " << E << " MeV" << G4 138 G4cout << "E " << E << " MeV" << G4endl;
141 G4cout << "mass " << E0 << " MeV" << G 139 G4cout << "mass " << E0 << " MeV" << G4endl;
142 G4cout << "charge " << Q << G4endl; 140 G4cout << "charge " << Q << G4endl;
143 } 141 }
144 Px = -Px; 142 Px = -Px;
145 Py = -Py; 143 Py = -Py;
146 Pz = -Pz; 144 Pz = -Pz;
147 145
148 G4double e1 = aParticle->GetKineticEnergy()/ 146 G4double e1 = aParticle->GetKineticEnergy()/MeV;
149 if (e1 < 1.) e1 = 1.; 147 if (e1 < 1.) e1 = 1.;
150 148
151 // Average number of neutrons 149 // Average number of neutrons
152 G4double avern = 2.569 + 0.559*G4Log(e1); 150 G4double avern = 2.569 + 0.559*G4Log(e1);
153 G4bool photofission = 0; // For now 151 G4bool photofission = 0; // For now
154 // Take the following value if photofission is 152 // Take the following value if photofission is not included
155 if (!photofission) avern = 2.569 + 0.900*G4 153 if (!photofission) avern = 2.569 + 0.900*G4Log(e1);
156 154
157 // Average number of gammas 155 // Average number of gammas
158 G4double averg = 9.500 + 0.600*G4Log(e1); 156 G4double averg = 9.500 + 0.600*G4Log(e1);
159 157
160 G4double ran = G4RandGauss::shoot(); 158 G4double ran = G4RandGauss::shoot();
161 // Number of neutrons 159 // Number of neutrons
162 G4int nn = static_cast<G4int>(avern + ran*1 160 G4int nn = static_cast<G4int>(avern + ran*1.23 + 0.5);
163 ran = G4RandGauss::shoot(); 161 ran = G4RandGauss::shoot();
164 // Number of gammas 162 // Number of gammas
165 G4int ng = static_cast<G4int>(averg + ran*3 163 G4int ng = static_cast<G4int>(averg + ran*3. + 0.5);
166 if (nn < 1) nn = 1; 164 if (nn < 1) nn = 1;
167 if (ng < 1) ng = 1; 165 if (ng < 1) ng = 1;
168 G4double exn = 0.; 166 G4double exn = 0.;
169 G4double exg = 0.; 167 G4double exg = 0.;
170 168
171 // Make secondary neutrons and distribute kine 169 // Make secondary neutrons and distribute kinetic energy
172 G4DynamicParticle* aNeutron; 170 G4DynamicParticle* aNeutron;
173 G4int i; 171 G4int i;
174 for (i = 1; i <= nn; i++) { 172 for (i = 1; i <= nn; i++) {
175 ran = G4UniformRand(); 173 ran = G4UniformRand();
176 G4int j; 174 G4int j;
177 for (j = 1; j <= 10; j++) { 175 for (j = 1; j <= 10; j++) {
178 if (ran < spneut[j-1]) goto label12; 176 if (ran < spneut[j-1]) goto label12;
179 } 177 }
180 j = 10; 178 j = 10;
181 label12: 179 label12:
182 ran = G4UniformRand(); 180 ran = G4UniformRand();
183 G4double ekin = (j - 1)*1. + ran; 181 G4double ekin = (j - 1)*1. + ran;
184 exn = exn + ekin; 182 exn = exn + ekin;
185 aNeutron = new G4DynamicParticle(G4Neutr 183 aNeutron = new G4DynamicParticle(G4Neutron::NeutronDefinition(),
186 G4Parti 184 G4ParticleMomentum(1.,0.,0.),
187 ekin*Me 185 ekin*MeV);
188 theParticleChange.AddSecondary(aNeutron, << 186 theParticleChange.AddSecondary(aNeutron);
189 } 187 }
190 188
191 // Make secondary gammas and distribute kineti 189 // Make secondary gammas and distribute kinetic energy
192 G4DynamicParticle* aGamma; 190 G4DynamicParticle* aGamma;
193 for (i = 1; i <= ng; i++) { 191 for (i = 1; i <= ng; i++) {
194 ran = G4UniformRand(); 192 ran = G4UniformRand();
195 G4double ekin = -0.87*G4Log(ran); 193 G4double ekin = -0.87*G4Log(ran);
196 exg = exg + ekin; 194 exg = exg + ekin;
197 aGamma = new G4DynamicParticle(G4Gamma:: 195 aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(),
198 G4Particl 196 G4ParticleMomentum(1.,0.,0.),
199 ekin*MeV) 197 ekin*MeV);
200 theParticleChange.AddSecondary(aGamma, s << 198 theParticleChange.AddSecondary(aGamma);
201 } 199 }
202 200
203 // Distribute momentum vectors and do Lorentz 201 // Distribute momentum vectors and do Lorentz transformation
204 202
205 G4HadSecondary* theSecondary; 203 G4HadSecondary* theSecondary;
206 204
207 for (i = 1; i <= nn + ng; i++) { 205 for (i = 1; i <= nn + ng; i++) {
208 G4double ran1 = G4UniformRand(); 206 G4double ran1 = G4UniformRand();
209 G4double ran2 = G4UniformRand(); 207 G4double ran2 = G4UniformRand();
210 G4double cost = -1. + 2.*ran1; 208 G4double cost = -1. + 2.*ran1;
211 G4double sint = std::sqrt(std::abs(1. - 209 G4double sint = std::sqrt(std::abs(1. - cost*cost));
212 G4double phi = ran2*twopi; 210 G4double phi = ran2*twopi;
213 // G4cout << ran1 << " " << ran2 << 211 // G4cout << ran1 << " " << ran2 << G4endl;
214 // G4cout << cost << " " << sint << 212 // G4cout << cost << " " << sint << " " << phi << G4endl;
215 theSecondary = theParticleChange.GetSeco 213 theSecondary = theParticleChange.GetSecondary(i - 1);
216 G4double pp = theSecondary->GetParticle( 214 G4double pp = theSecondary->GetParticle()->GetTotalMomentum()/MeV;
217 G4double px = pp*sint*std::sin(phi); 215 G4double px = pp*sint*std::sin(phi);
218 G4double py = pp*sint*std::cos(phi); 216 G4double py = pp*sint*std::cos(phi);
219 G4double pz = pp*cost; 217 G4double pz = pp*cost;
220 // G4cout << pp << G4endl; 218 // G4cout << pp << G4endl;
221 // G4cout << px << " " << py << " " 219 // G4cout << px << " " << py << " " << pz << G4endl;
222 G4double e = theSecondary->GetParticle() 220 G4double e = theSecondary->GetParticle()->GetTotalEnergy()/MeV;
223 G4double e0 = theSecondary->GetParticle( 221 G4double e0 = theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV;
224 222
225 G4double a = px*Px + py*Py + pz*Pz; 223 G4double a = px*Px + py*Py + pz*Pz;
226 a = (a/(E + E0) - e)/E0; 224 a = (a/(E + E0) - e)/E0;
227 225
228 px = px + a*Px; 226 px = px + a*Px;
229 py = py + a*Py; 227 py = py + a*Py;
230 pz = pz + a*Pz; 228 pz = pz + a*Pz;
231 G4double p2 = px*px + py*py + pz*pz; 229 G4double p2 = px*px + py*py + pz*pz;
232 pp = std::sqrt(p2); 230 pp = std::sqrt(p2);
233 e = std::sqrt(e0*e0 + p2); 231 e = std::sqrt(e0*e0 + p2);
234 G4double ekin = e - theSecondary->GetPar 232 G4double ekin = e - theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV;
235 theSecondary->GetParticle()->SetMomentum 233 theSecondary->GetParticle()->SetMomentumDirection(G4ParticleMomentum(px/pp,
236 234 py/pp,
237 235 pz/pp));
238 theSecondary->GetParticle()->SetKineticE 236 theSecondary->GetParticle()->SetKineticEnergy(ekin*MeV);
239 } 237 }
240 238
241 return &theParticleChange; 239 return &theParticleChange;
242 } 240 }
243 241
244 // Computes atomic mass in MeV (translation of 242 // Computes atomic mass in MeV (translation of GHEISHA routine ATOMAS)
245 // Not optimized: conforms closely to original 243 // Not optimized: conforms closely to original Fortran.
246 244
247 G4double G4LFission::Atomas(const G4double A, 245 G4double G4LFission::Atomas(const G4double A, const G4double Z)
248 { 246 {
249 G4double rmel = G4Electron::ElectronDefiniti 247 G4double rmel = G4Electron::ElectronDefinition()->GetPDGMass()/MeV;
250 G4double rmp = G4Proton::ProtonDefinition() 248 G4double rmp = G4Proton::ProtonDefinition()->GetPDGMass()/MeV;
251 G4double rmn = G4Neutron::NeutronDefinition 249 G4double rmn = G4Neutron::NeutronDefinition()->GetPDGMass()/MeV;
252 G4double rmd = G4Deuteron::DeuteronDefiniti 250 G4double rmd = G4Deuteron::DeuteronDefinition()->GetPDGMass()/MeV;
253 G4double rma = G4Alpha::AlphaDefinition()-> 251 G4double rma = G4Alpha::AlphaDefinition()->GetPDGMass()/MeV;
254 252
255 G4int ia = static_cast<G4int>(A + 0.5); 253 G4int ia = static_cast<G4int>(A + 0.5);
256 if (ia < 1) return 0; 254 if (ia < 1) return 0;
257 G4int iz = static_cast<G4int>(Z + 0.5); 255 G4int iz = static_cast<G4int>(Z + 0.5);
258 if (iz < 0) return 0; 256 if (iz < 0) return 0;
259 if (iz > ia) return 0; 257 if (iz > ia) return 0;
260 258
261 if (ia == 1) { 259 if (ia == 1) {
262 if (iz == 0) return rmn; //neut 260 if (iz == 0) return rmn; //neutron
263 if (iz == 1) return rmp + rmel; //Hydr 261 if (iz == 1) return rmp + rmel; //Hydrogen
264 } 262 }
265 else if (ia == 2 && iz == 1) { 263 else if (ia == 2 && iz == 1) {
266 return rmd; //Deut 264 return rmd; //Deuteron
267 } 265 }
268 else if (ia == 4 && iz == 2) { 266 else if (ia == 4 && iz == 2) {
269 return rma; //Alph 267 return rma; //Alpha
270 } 268 }
271 269
272 G4Pow* Pow=G4Pow::GetInstance(); 270 G4Pow* Pow=G4Pow::GetInstance();
273 G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel 271 G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel - 15.67*A
274 + 17.23*Pow->A23(A) 272 + 17.23*Pow->A23(A)
275 + 93.15*(A/2. - Z)*(A/2. - Z 273 + 93.15*(A/2. - Z)*(A/2. - Z)/A
276 + 0.6984523*Z*Z/Pow->A13(A); 274 + 0.6984523*Z*Z/Pow->A13(A);
277 G4int ipp = (ia - iz)%2; 275 G4int ipp = (ia - iz)%2;
278 G4int izz = iz%2; 276 G4int izz = iz%2;
279 if (ipp == izz) mass = mass + (ipp + izz -1) 277 if (ipp == izz) mass = mass + (ipp + izz -1)*12.*Pow->powA(A, -0.5);
280 278
281 return mass; 279 return mass;
282 } 280 }
283 281
284 const std::pair<G4double, G4double> G4LFission 282 const std::pair<G4double, G4double> G4LFission::GetFatalEnergyCheckLevels() const
285 { 283 {
286 // max energy non-conservation is mass of he 284 // max energy non-conservation is mass of heavy nucleus
287 return std::pair<G4double, G4double>(10.0*pe << 285 return std::pair<G4double, G4double>(5*perCent,250*GeV);
288 } 286 }
289 287