Geant4  10.00.p03
G4LFission.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 //
26 // $Id: G4LFission.cc 66892 2013-01-17 10:57:59Z gunter $
27 //
28 //
29 // G4 Model: Low Energy Fission
30 // F.W. Jones, TRIUMF, 03-DEC-96
31 //
32 // This is a prototype of a low-energy fission process.
33 // Currently it is based on the GHEISHA routine FISSIO,
34 // and conforms fairly closely to the original Fortran.
35 // Note: energy is in MeV and momentum is in MeV/c.
36 //
37 // use -scheme for elastic scattering: HPW, 20th June 1997
38 // the code comes mostly from the old Low-energy Fission class
39 //
40 // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
41 
42 #include <iostream>
43 
44 #include "G4LFission.hh"
45 #include "globals.hh"
46 #include "G4PhysicalConstants.hh"
47 #include "G4SystemOfUnits.hh"
48 #include "Randomize.hh"
49 
51  : G4HadronicInteraction(name)
52 {
53  init();
54  SetMinEnergy(0.0*GeV);
56 }
57 
58 
60 {
62 }
63 
64 
65 void G4LFission::ModelDescription(std::ostream& outFile) const
66 {
67  outFile << "G4LFission is one of the Low Energy Parameterized\n"
68  << "(LEP) models used to implement neutron-induced fission of\n"
69  << "nuclei. It is a re-engineered version of the GHEISHA code\n"
70  << "of H. Fesefeldt which emits neutrons and gammas but no\n"
71  << "nuclear fragments. The model is applicable to all incident\n"
72  << "neutron energies.\n";
73 }
74 
76 {
77  G4int i;
78  G4double xx = 1. - 0.5;
79  G4double xxx = std::sqrt(2.29*xx);
80  spneut[0] = std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.;
81  for (i = 2; i <= 10; i++) {
82  xx = i*1. - 0.5;
83  xxx = std::sqrt(2.29*xx);
84  spneut[i-1] = spneut[i-2] + std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.;
85  }
86  for (i = 1; i <= 10; i++) {
87  spneut[i-1] = spneut[i-1]/spneut[9];
88  if (verboseLevel > 1) G4cout << "G4LFission::init: i=" << i <<
89  " spneut=" << spneut[i-1] << G4endl;
90  }
91 }
92 
93 
95  G4Nucleus& targetNucleus)
96 {
98  const G4HadProjectile* aParticle = &aTrack;
99 
100  G4double N = targetNucleus.GetA_asInt();
101  G4double Z = targetNucleus.GetZ_asInt();
103 
104  G4double P = aParticle->GetTotalMomentum()/MeV;
105  G4double Px = aParticle->Get4Momentum().vect().x();
106  G4double Py = aParticle->Get4Momentum().vect().y();
107  G4double Pz = aParticle->Get4Momentum().vect().z();
108  G4double E = aParticle->GetTotalEnergy()/MeV;
109  G4double E0 = aParticle->GetDefinition()->GetPDGMass()/MeV;
110  G4double Q = aParticle->GetDefinition()->GetPDGCharge();
111  if (verboseLevel > 1) {
112  G4cout << "G4LFission:ApplyYourself: incident particle:" << G4endl;
113  G4cout << "P " << P << " MeV/c" << G4endl;
114  G4cout << "Px " << Px << " MeV/c" << G4endl;
115  G4cout << "Py " << Py << " MeV/c" << G4endl;
116  G4cout << "Pz " << Pz << " MeV/c" << G4endl;
117  G4cout << "E " << E << " MeV" << G4endl;
118  G4cout << "mass " << E0 << " MeV" << G4endl;
119  G4cout << "charge " << Q << G4endl;
120  }
121  // GHEISHA ADD operation to get total energy, mass, charge:
122  if (verboseLevel > 1) {
123  G4cout << "G4LFission:ApplyYourself: material:" << G4endl;
124  G4cout << "A " << N << G4endl;
125  G4cout << "Z " << Z << G4endl;
126  G4cout << "atomic mass " <<
127  Atomas(N, Z) << "MeV" << G4endl;
128  }
129  E = E + Atomas(N, Z);
130  G4double E02 = E*E - P*P;
131  E0 = std::sqrt(std::abs(E02));
132  if (E02 < 0) E0 = -E0;
133  Q = Q + Z;
134  if (verboseLevel > 1) {
135  G4cout << "G4LFission:ApplyYourself: total:" << G4endl;
136  G4cout << "E " << E << " MeV" << G4endl;
137  G4cout << "mass " << E0 << " MeV" << G4endl;
138  G4cout << "charge " << Q << G4endl;
139  }
140  Px = -Px;
141  Py = -Py;
142  Pz = -Pz;
143 
144  G4double e1 = aParticle->GetKineticEnergy()/MeV;
145  if (e1 < 1.) e1 = 1.;
146 
147 // Average number of neutrons
148  G4double avern = 2.569 + 0.559*std::log(e1);
149  G4bool photofission = 0; // For now
150 // Take the following value if photofission is not included
151  if (!photofission) avern = 2.569 + 0.900*std::log(e1);
152 
153 // Average number of gammas
154  G4double averg = 9.500 + 0.600*std::log(e1);
155 
157 // Number of neutrons
158  G4int nn = static_cast<G4int>(avern + ran*1.23 + 0.5);
159  ran = G4RandGauss::shoot();
160 // Number of gammas
161  G4int ng = static_cast<G4int>(averg + ran*3. + 0.5);
162  if (nn < 1) nn = 1;
163  if (ng < 1) ng = 1;
164  G4double exn = 0.;
165  G4double exg = 0.;
166 
167 // Make secondary neutrons and distribute kinetic energy
168  G4DynamicParticle* aNeutron;
169  G4int i;
170  for (i = 1; i <= nn; i++) {
171  ran = G4UniformRand();
172  G4int j;
173  for (j = 1; j <= 10; j++) {
174  if (ran < spneut[j-1]) goto label12;
175  }
176  j = 10;
177  label12:
178  ran = G4UniformRand();
179  G4double ekin = (j - 1)*1. + ran;
180  exn = exn + ekin;
182  G4ParticleMomentum(1.,0.,0.),
183  ekin*MeV);
185  }
186 
187 // Make secondary gammas and distribute kinetic energy
188  G4DynamicParticle* aGamma;
189  for (i = 1; i <= ng; i++) {
190  ran = G4UniformRand();
191  G4double ekin = -0.87*std::log(ran);
192  exg = exg + ekin;
194  G4ParticleMomentum(1.,0.,0.),
195  ekin*MeV);
197  }
198 
199 // Distribute momentum vectors and do Lorentz transformation
200 
201  G4HadSecondary* theSecondary;
202 
203  for (i = 1; i <= nn + ng; i++) {
204  G4double ran1 = G4UniformRand();
205  G4double ran2 = G4UniformRand();
206  G4double cost = -1. + 2.*ran1;
207  G4double sint = std::sqrt(std::abs(1. - cost*cost));
208  G4double phi = ran2*twopi;
209  // G4cout << ran1 << " " << ran2 << G4endl;
210  // G4cout << cost << " " << sint << " " << phi << G4endl;
211  theSecondary = theParticleChange.GetSecondary(i - 1);
212  G4double pp = theSecondary->GetParticle()->GetTotalMomentum()/MeV;
213  G4double px = pp*sint*std::sin(phi);
214  G4double py = pp*sint*std::cos(phi);
215  G4double pz = pp*cost;
216  // G4cout << pp << G4endl;
217  // G4cout << px << " " << py << " " << pz << G4endl;
218  G4double e = theSecondary->GetParticle()->GetTotalEnergy()/MeV;
219  G4double e0 = theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV;
220 
221  G4double a = px*Px + py*Py + pz*Pz;
222  a = (a/(E + E0) - e)/E0;
223 
224  px = px + a*Px;
225  py = py + a*Py;
226  pz = pz + a*Pz;
227  G4double p2 = px*px + py*py + pz*pz;
228  pp = std::sqrt(p2);
229  e = std::sqrt(e0*e0 + p2);
230  G4double ekin = e - theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV;
231  theSecondary->GetParticle()->SetMomentumDirection(G4ParticleMomentum(px/pp,
232  py/pp,
233  pz/pp));
234  theSecondary->GetParticle()->SetKineticEnergy(ekin*MeV);
235  }
236 
237  return &theParticleChange;
238 }
239 
240 // Computes atomic mass in MeV (translation of GHEISHA routine ATOMAS)
241 // Not optimized: conforms closely to original Fortran.
242 
244 {
250 
251  G4int ia = static_cast<G4int>(A + 0.5);
252  if (ia < 1) return 0;
253  G4int iz = static_cast<G4int>(Z + 0.5);
254  if (iz < 0) return 0;
255  if (iz > ia) return 0;
256 
257  if (ia == 1) {
258  if (iz == 0) return rmn; //neutron
259  if (iz == 1) return rmp + rmel; //Hydrogen
260  }
261  else if (ia == 2 && iz == 1) {
262  return rmd; //Deuteron
263  }
264  else if (ia == 4 && iz == 2) {
265  return rma; //Alpha
266  }
267 
268  G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel - 15.67*A
269  + 17.23*std::pow(A, 2./3.)
270  + 93.15*(A/2. - Z)*(A/2. - Z)/A
271  + 0.6984523*Z*Z/std::pow(A, 1./3.);
272  G4int ipp = (ia - iz)%2;
273  G4int izz = iz%2;
274  if (ipp == izz) mass = mass + (ipp + izz -1)*12.*std::pow(A, -0.5);
275 
276  return mass;
277 }
278 
279 const std::pair<G4double, G4double> G4LFission::GetFatalEnergyCheckLevels() const
280 {
281  // max energy non-conservation is mass of heavy nucleus
282  return std::pair<G4double, G4double>(5*perCent,250*GeV);
283 }
static G4Electron * ElectronDefinition()
Definition: G4Electron.cc:89
G4int GetA_asInt() const
Definition: G4Nucleus.hh:109
ThreeVector shoot(const G4int Ap, const G4int Af)
static const double MeV
Definition: G4SIunits.hh:193
G4HadSecondary * GetSecondary(size_t i)
virtual const std::pair< G4double, G4double > GetFatalEnergyCheckLevels() const
Definition: G4LFission.cc:279
G4double GetTotalEnergy() const
std::ofstream outFile
Definition: GammaRayTel.cc:68
virtual void ModelDescription(std::ostream &outFile) const
Definition: G4LFission.cc:65
G4String name
Definition: TRTMaterials.hh:40
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4ParticleDefinition * GetDefinition() const
G4double a
Definition: TRTMaterials.hh:39
static G4double Atomas(const G4double A, const G4double Z)
Definition: G4LFission.cc:243
int G4int
Definition: G4Types.hh:78
G4double GetTotalMomentum() const
void SetStatusChange(G4HadFinalStateStatus aS)
void SetMinEnergy(G4double anEnergy)
void init()
Definition: G4LFission.cc:75
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
const G4ParticleDefinition * GetDefinition() const
bool G4bool
Definition: G4Types.hh:79
G4double iz
Definition: TRTMaterials.hh:39
G4double GetKineticEnergy() const
G4LFission(const G4String &name="G4LFission")
Definition: G4LFission.cc:50
static const double GeV
Definition: G4SIunits.hh:196
static const double perCent
Definition: G4SIunits.hh:296
static const G4double A[nN]
const G4LorentzVector & Get4Momentum() const
G4double spneut[10]
Definition: G4LFission.hh:84
void SetKineticEnergy(G4double aEnergy)
static const G4double e1
G4HadFinalState * ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
Definition: G4LFission.cc:94
G4double GetPDGMass() const
G4DynamicParticle * GetParticle()
G4int GetZ_asInt() const
Definition: G4Nucleus.hh:115
void SetMaxEnergy(const G4double anEnergy)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
G4double GetPDGCharge() const
G4ThreeVector G4ParticleMomentum
#define DBL_MAX
Definition: templates.hh:83
static G4Deuteron * DeuteronDefinition()
Definition: G4Deuteron.cc:89
static G4Alpha * AlphaDefinition()
Definition: G4Alpha.cc:84
static G4Neutron * NeutronDefinition()
Definition: G4Neutron.cc:99
void AddSecondary(G4DynamicParticle *aP)
G4double GetTotalMomentum() const
G4double GetTotalEnergy() const
static G4Gamma * GammaDefinition()
Definition: G4Gamma.cc:81