// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // $Id: G4LFission.cc,v 1.15 2007/02/26 19:29:30 dennis Exp $ // GEANT4 tag $Name: geant4-09-02 $ // // // G4 Model: Low Energy Fission // F.W. Jones, TRIUMF, 03-DEC-96 // // This is a prototype of a low-energy fission process. // Currently it is based on the GHEISHA routine FISSIO, // and conforms fairly closely to the original Fortran. // Note: energy is in MeV and momentum is in MeV/c. // // use -scheme for elastic scattering: HPW, 20th June 1997 // the code comes mostly from the old Low-energy Fission class // // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange. // #include "globals.hh" #include "G4LFission.hh" #include "Randomize.hh" G4LFission::G4LFission() : G4HadronicInteraction("G4LFission") { init(); SetMinEnergy( 0.0*GeV ); SetMaxEnergy( DBL_MAX ); } G4LFission::~G4LFission() { theParticleChange.Clear(); } void G4LFission::init() { G4int i; G4double xx = 1. - 0.5; G4double xxx = std::sqrt(2.29*xx); spneut[0] = std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.; for (i = 2; i <= 10; i++) { xx = i*1. - 0.5; xxx = std::sqrt(2.29*xx); spneut[i-1] = spneut[i-2] + std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.; } for (i = 1; i <= 10; i++) { spneut[i-1] = spneut[i-1]/spneut[9]; if (verboseLevel > 1) G4cout << "G4LFission::init: i=" << i << " spneut=" << spneut[i-1] << G4endl; } } G4HadFinalState* G4LFission::ApplyYourself(const G4HadProjectile & aTrack,G4Nucleus & targetNucleus) { theParticleChange.Clear(); const G4HadProjectile* aParticle = &aTrack; G4double N = targetNucleus.GetN(); G4double Z = targetNucleus.GetZ(); theParticleChange.SetStatusChange(stopAndKill); G4double P = aParticle->GetTotalMomentum()/MeV; G4double Px = aParticle->Get4Momentum().vect().x(); G4double Py = aParticle->Get4Momentum().vect().y(); G4double Pz = aParticle->Get4Momentum().vect().z(); G4double E = aParticle->GetTotalEnergy()/MeV; G4double E0 = aParticle->GetDefinition()->GetPDGMass()/MeV; G4double Q = aParticle->GetDefinition()->GetPDGCharge(); if (verboseLevel > 1) { G4cout << "G4LFission:ApplyYourself: incident particle:" << G4endl; G4cout << "P " << P << " MeV/c" << G4endl; G4cout << "Px " << Px << " MeV/c" << G4endl; G4cout << "Py " << Py << " MeV/c" << G4endl; G4cout << "Pz " << Pz << " MeV/c" << G4endl; G4cout << "E " << E << " MeV" << G4endl; G4cout << "mass " << E0 << " MeV" << G4endl; G4cout << "charge " << Q << G4endl; } // GHEISHA ADD operation to get total energy, mass, charge: if (verboseLevel > 1) { G4cout << "G4LFission:ApplyYourself: material:" << G4endl; G4cout << "A " << N << G4endl; G4cout << "Z " << Z << G4endl; G4cout << "atomic mass " << Atomas(N, Z) << "MeV" << G4endl; } E = E + Atomas(N, Z); G4double E02 = E*E - P*P; E0 = std::sqrt(std::abs(E02)); if (E02 < 0) E0 = -E0; Q = Q + Z; if (verboseLevel > 1) { G4cout << "G4LFission:ApplyYourself: total:" << G4endl; G4cout << "E " << E << " MeV" << G4endl; G4cout << "mass " << E0 << " MeV" << G4endl; G4cout << "charge " << Q << G4endl; } Px = -Px; Py = -Py; Pz = -Pz; G4double e1 = aParticle->GetKineticEnergy()/MeV; if (e1 < 1.) e1 = 1.; // Average number of neutrons G4double avern = 2.569 + 0.559*std::log(e1); G4bool photofission = 0; // For now // Take the following value if photofission is not included if (!photofission) avern = 2.569 + 0.900*std::log(e1); // Average number of gammas G4double averg = 9.500 + 0.600*std::log(e1); G4double ran = G4RandGauss::shoot(); // Number of neutrons G4int nn = static_cast(avern + ran*1.23 + 0.5); ran = G4RandGauss::shoot(); // Number of gammas G4int ng = static_cast(averg + ran*3. + 0.5); if (nn < 1) nn = 1; if (ng < 1) ng = 1; G4double exn = 0.; G4double exg = 0.; // Make secondary neutrons and distribute kinetic energy G4DynamicParticle* aNeutron; G4int i; for (i = 1; i <= nn; i++) { ran = G4UniformRand(); G4int j; for (j = 1; j <= 10; j++) { if (ran < spneut[j-1]) goto label12; } j = 10; label12: ran = G4UniformRand(); G4double ekin = (j - 1)*1. + ran; exn = exn + ekin; aNeutron = new G4DynamicParticle(G4Neutron::NeutronDefinition(), G4ParticleMomentum(1.,0.,0.), ekin*MeV); theParticleChange.AddSecondary(aNeutron); } // Make secondary gammas and distribute kinetic energy G4DynamicParticle* aGamma; for (i = 1; i <= ng; i++) { ran = G4UniformRand(); G4double ekin = -0.87*std::log(ran); exg = exg + ekin; aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(), G4ParticleMomentum(1.,0.,0.), ekin*MeV); theParticleChange.AddSecondary(aGamma); } // Distribute momentum vectors and do Lorentz transformation G4HadSecondary* theSecondary; for (i = 1; i <= nn + ng; i++) { G4double ran1 = G4UniformRand(); G4double ran2 = G4UniformRand(); G4double cost = -1. + 2.*ran1; G4double sint = std::sqrt(std::abs(1. - cost*cost)); G4double phi = ran2*twopi; // G4cout << ran1 << " " << ran2 << G4endl; // G4cout << cost << " " << sint << " " << phi << G4endl; theSecondary = theParticleChange.GetSecondary(i - 1); G4double pp = theSecondary->GetParticle()->GetTotalMomentum()/MeV; G4double px = pp*sint*std::sin(phi); G4double py = pp*sint*std::cos(phi); G4double pz = pp*cost; // G4cout << pp << G4endl; // G4cout << px << " " << py << " " << pz << G4endl; G4double e = theSecondary->GetParticle()->GetTotalEnergy()/MeV; G4double e0 = theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV; G4double a = px*Px + py*Py + pz*Pz; a = (a/(E + E0) - e)/E0; px = px + a*Px; py = py + a*Py; pz = pz + a*Pz; G4double p2 = px*px + py*py + pz*pz; pp = std::sqrt(p2); e = std::sqrt(e0*e0 + p2); G4double ekin = e - theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV; theSecondary->GetParticle()->SetMomentumDirection(G4ParticleMomentum(px/pp, py/pp, pz/pp)); theSecondary->GetParticle()->SetKineticEnergy(ekin*MeV); } return &theParticleChange; } // Computes atomic mass in MeV (translation of GHEISHA routine ATOMAS) // Not optimized: conforms closely to original Fortran. G4double G4LFission::Atomas(const G4double A, const G4double Z) { G4double rmel = G4Electron::ElectronDefinition()->GetPDGMass()/MeV; G4double rmp = G4Proton::ProtonDefinition()->GetPDGMass()/MeV; G4double rmn = G4Neutron::NeutronDefinition()->GetPDGMass()/MeV; G4double rmd = G4Deuteron::DeuteronDefinition()->GetPDGMass()/MeV; G4double rma = G4Alpha::AlphaDefinition()->GetPDGMass()/MeV; G4int ia = static_cast(A + 0.5); if (ia < 1) return 0; G4int iz = static_cast(Z + 0.5); if (iz < 0) return 0; if (iz > ia) return 0; if (ia == 1) { if (iz == 0) return rmn; //neutron if (iz == 1) return rmp + rmel; //Hydrogen } else if (ia == 2 && iz == 1) { return rmd; //Deuteron } else if (ia == 4 && iz == 2) { return rma; //Alpha } G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel - 15.67*A + 17.23*std::pow(A, 2./3.) + 93.15*(A/2. - Z)*(A/2. - Z)/A + 0.6984523*Z*Z/std::pow(A, 1./3.); G4int ipp = (ia - iz)%2; G4int izz = iz%2; if (ipp == izz) mass = mass + (ipp + izz -1)*12.*std::pow(A, -0.5); return mass; }