source: trunk/source/processes/hadronic/models/low_energy/src/G4LFission.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 9.5 KB
Line 
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-03-cand-01 $
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
49G4LFission::G4LFission() : G4HadronicInteraction("G4LFission")
50{
51 init();
52 SetMinEnergy( 0.0*GeV );
53 SetMaxEnergy( DBL_MAX );
54
55}
56
57G4LFission::~G4LFission()
58{
59 theParticleChange.Clear();
60}
61
62
63void
64G4LFission::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
83G4HadFinalState*
84G4LFission::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
232G4double
233G4LFission::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}
Note: See TracBrowser for help on using the repository browser.