G4LFission.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 //
00029 // G4 Model: Low Energy Fission
00030 // F.W. Jones, TRIUMF, 03-DEC-96
00031 // 
00032 // This is a prototype of a low-energy fission process.
00033 // Currently it is based on the GHEISHA routine FISSIO,
00034 // and conforms fairly closely to the original Fortran.
00035 // Note: energy is in MeV and momentum is in MeV/c.
00036 //
00037 // use -scheme for elastic scattering: HPW, 20th June 1997
00038 // the code comes mostly from the old Low-energy Fission class
00039 //
00040 // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
00041 
00042 #include <iostream>
00043 
00044 #include "G4LFission.hh"
00045 #include "globals.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "Randomize.hh"
00049 
00050 G4LFission::G4LFission(const G4String& name)
00051  : G4HadronicInteraction(name)
00052 {
00053   init();
00054   SetMinEnergy(0.0*GeV);
00055   SetMaxEnergy(DBL_MAX);
00056 }
00057 
00058 
00059 G4LFission::~G4LFission()
00060 {
00061   theParticleChange.Clear();
00062 }
00063 
00064 
00065 void G4LFission::ModelDescription(std::ostream& outFile) const
00066 {
00067   outFile << "G4LFission is one of the Low Energy Parameterized\n"
00068           << "(LEP) models used to implement neutron-induced fission of\n"
00069           << "nuclei.  It is a re-engineered version of the GHEISHA code\n"
00070           << "of H. Fesefeldt which emits neutrons and gammas but no\n"
00071           << "nuclear fragments.  The model is applicable to all incident\n"
00072           << "neutron energies.\n";
00073 }
00074 
00075 void G4LFission::init()
00076 {
00077    G4int i;
00078    G4double xx = 1. - 0.5;
00079    G4double xxx = std::sqrt(2.29*xx);
00080    spneut[0] = std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.;
00081    for (i = 2; i <= 10; i++) {
00082       xx = i*1. - 0.5;
00083       xxx = std::sqrt(2.29*xx);
00084       spneut[i-1] = spneut[i-2] + std::exp(-xx/0.965)*(std::exp(xxx) - std::exp(-xxx))/2.;
00085    }
00086    for (i = 1; i <= 10; i++) {
00087       spneut[i-1] = spneut[i-1]/spneut[9];
00088       if (verboseLevel > 1) G4cout << "G4LFission::init: i=" << i << 
00089          " spneut=" << spneut[i-1] << G4endl;
00090    }
00091 }
00092 
00093 
00094 G4HadFinalState* G4LFission::ApplyYourself(const G4HadProjectile& aTrack,
00095                                            G4Nucleus& targetNucleus)
00096 {
00097   theParticleChange.Clear();
00098   const G4HadProjectile* aParticle = &aTrack;
00099 
00100   G4double N = targetNucleus.GetA_asInt();
00101   G4double Z = targetNucleus.GetZ_asInt();
00102   theParticleChange.SetStatusChange(stopAndKill);
00103 
00104   G4double P = aParticle->GetTotalMomentum()/MeV;
00105   G4double Px = aParticle->Get4Momentum().vect().x();
00106   G4double Py = aParticle->Get4Momentum().vect().y();
00107   G4double Pz = aParticle->Get4Momentum().vect().z();
00108   G4double E = aParticle->GetTotalEnergy()/MeV;
00109   G4double E0 = aParticle->GetDefinition()->GetPDGMass()/MeV;
00110   G4double Q = aParticle->GetDefinition()->GetPDGCharge();
00111   if (verboseLevel > 1) {
00112       G4cout << "G4LFission:ApplyYourself: incident particle:" << G4endl;
00113       G4cout << "P      " << P << " MeV/c" << G4endl;
00114       G4cout << "Px     " << Px << " MeV/c" << G4endl;
00115       G4cout << "Py     " << Py << " MeV/c" << G4endl;
00116       G4cout << "Pz     " << Pz << " MeV/c" << G4endl;
00117       G4cout << "E      " << E << " MeV" << G4endl;
00118       G4cout << "mass   " << E0 << " MeV" << G4endl;
00119       G4cout << "charge " << Q << G4endl;
00120   }
00121   // GHEISHA ADD operation to get total energy, mass, charge:
00122    if (verboseLevel > 1) {
00123       G4cout << "G4LFission:ApplyYourself: material:" << G4endl;
00124       G4cout << "A      " << N << G4endl;
00125       G4cout << "Z      " << Z << G4endl;
00126       G4cout << "atomic mass " << 
00127         Atomas(N, Z) << "MeV" << G4endl;
00128    }
00129   E = E + Atomas(N, Z);
00130   G4double E02 = E*E - P*P;
00131   E0 = std::sqrt(std::abs(E02));
00132   if (E02 < 0) E0 = -E0;
00133   Q = Q + Z;
00134   if (verboseLevel > 1) {
00135       G4cout << "G4LFission:ApplyYourself: total:" << G4endl;
00136       G4cout << "E      " << E << " MeV" << G4endl;
00137       G4cout << "mass   " << E0 << " MeV" << G4endl;
00138       G4cout << "charge " << Q << G4endl;
00139   }
00140   Px = -Px;
00141   Py = -Py;
00142   Pz = -Pz;
00143 
00144   G4double e1 = aParticle->GetKineticEnergy()/MeV;
00145    if (e1 < 1.) e1 = 1.;
00146 
00147 // Average number of neutrons
00148    G4double avern = 2.569 + 0.559*std::log(e1);
00149    G4bool photofission = 0;      // For now
00150 // Take the following value if photofission is not included
00151    if (!photofission) avern = 2.569 + 0.900*std::log(e1);
00152 
00153 // Average number of gammas
00154    G4double averg = 9.500 + 0.600*std::log(e1);
00155 
00156    G4double ran = G4RandGauss::shoot();
00157 // Number of neutrons
00158    G4int nn = static_cast<G4int>(avern + ran*1.23 + 0.5);
00159    ran = G4RandGauss::shoot();
00160 // Number of gammas
00161    G4int ng = static_cast<G4int>(averg + ran*3. + 0.5);
00162    if (nn < 1) nn = 1;
00163    if (ng < 1) ng = 1;
00164    G4double exn = 0.;
00165    G4double exg = 0.;
00166 
00167 // Make secondary neutrons and distribute kinetic energy
00168    G4DynamicParticle* aNeutron;
00169    G4int i;
00170    for (i = 1; i <= nn; i++) {
00171       ran = G4UniformRand();
00172       G4int j;
00173       for (j = 1; j <= 10; j++) {
00174          if (ran < spneut[j-1]) goto label12;
00175       }
00176       j = 10;
00177     label12:
00178       ran = G4UniformRand();
00179       G4double ekin = (j - 1)*1. + ran;
00180       exn = exn + ekin;
00181       aNeutron = new G4DynamicParticle(G4Neutron::NeutronDefinition(),
00182                                        G4ParticleMomentum(1.,0.,0.),
00183                                        ekin*MeV);
00184       theParticleChange.AddSecondary(aNeutron);
00185    }
00186 
00187 // Make secondary gammas and distribute kinetic energy
00188    G4DynamicParticle* aGamma;
00189    for (i = 1; i <= ng; i++) {
00190       ran = G4UniformRand();
00191       G4double ekin = -0.87*std::log(ran);
00192       exg = exg + ekin;
00193       aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(),
00194                                      G4ParticleMomentum(1.,0.,0.),
00195                                      ekin*MeV);
00196       theParticleChange.AddSecondary(aGamma);
00197    }
00198 
00199 // Distribute momentum vectors and do Lorentz transformation
00200 
00201    G4HadSecondary* theSecondary;
00202 
00203    for (i = 1; i <= nn + ng; i++) {
00204       G4double ran1 = G4UniformRand();
00205       G4double ran2 = G4UniformRand();
00206       G4double cost = -1. + 2.*ran1;
00207       G4double sint = std::sqrt(std::abs(1. - cost*cost));
00208       G4double phi = ran2*twopi;
00209       //      G4cout << ran1 << " " << ran2 << G4endl;
00210       //      G4cout << cost << " " << sint << " " << phi << G4endl;
00211       theSecondary = theParticleChange.GetSecondary(i - 1);
00212       G4double pp = theSecondary->GetParticle()->GetTotalMomentum()/MeV;
00213       G4double px = pp*sint*std::sin(phi);
00214       G4double py = pp*sint*std::cos(phi);
00215       G4double pz = pp*cost;
00216       //      G4cout << pp << G4endl;
00217       //      G4cout << px << " " << py << " " << pz << G4endl;
00218       G4double e = theSecondary->GetParticle()->GetTotalEnergy()/MeV;
00219       G4double e0 = theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV;
00220 
00221       G4double a = px*Px + py*Py + pz*Pz;
00222       a = (a/(E + E0) - e)/E0;
00223 
00224       px = px + a*Px;
00225       py = py + a*Py;
00226       pz = pz + a*Pz;
00227       G4double p2 = px*px + py*py + pz*pz;
00228       pp = std::sqrt(p2);
00229       e = std::sqrt(e0*e0 + p2);
00230       G4double ekin = e - theSecondary->GetParticle()->GetDefinition()->GetPDGMass()/MeV;
00231       theSecondary->GetParticle()->SetMomentumDirection(G4ParticleMomentum(px/pp,
00232                                                             py/pp,
00233                                                             pz/pp));
00234       theSecondary->GetParticle()->SetKineticEnergy(ekin*MeV);
00235    }
00236    
00237   return &theParticleChange;
00238 }
00239 
00240 // Computes atomic mass in MeV (translation of GHEISHA routine ATOMAS)
00241 // Not optimized: conforms closely to original Fortran.
00242 
00243 G4double G4LFission::Atomas(const G4double A, const G4double Z)
00244 {
00245   G4double rmel = G4Electron::ElectronDefinition()->GetPDGMass()/MeV;
00246   G4double rmp  = G4Proton::ProtonDefinition()->GetPDGMass()/MeV;
00247   G4double rmn  = G4Neutron::NeutronDefinition()->GetPDGMass()/MeV;
00248   G4double rmd  = G4Deuteron::DeuteronDefinition()->GetPDGMass()/MeV;
00249   G4double rma  = G4Alpha::AlphaDefinition()->GetPDGMass()/MeV;
00250 
00251   G4int ia = static_cast<G4int>(A + 0.5);
00252    if (ia < 1) return 0;
00253    G4int iz = static_cast<G4int>(Z + 0.5);
00254    if (iz < 0) return 0;
00255    if (iz > ia) return 0;
00256 
00257    if (ia == 1) {
00258       if (iz == 0) return rmn;          //neutron
00259       if (iz == 1) return rmp + rmel;   //Hydrogen
00260    }
00261    else if (ia == 2 && iz == 1) {
00262       return rmd;                       //Deuteron
00263    }
00264    else if (ia == 4 && iz == 2) {
00265       return rma;                       //Alpha
00266    }
00267 
00268   G4double mass = (A - Z)*rmn + Z*rmp + Z*rmel - 15.67*A
00269                   + 17.23*std::pow(A, 2./3.)
00270                   + 93.15*(A/2. - Z)*(A/2. - Z)/A
00271                   + 0.6984523*Z*Z/std::pow(A, 1./3.);
00272   G4int ipp = (ia - iz)%2;
00273   G4int izz = iz%2;
00274   if (ipp == izz) mass = mass + (ipp + izz -1)*12.*std::pow(A, -0.5);
00275 
00276   return mass;
00277 }
00278 
00279 const std::pair<G4double, G4double> G4LFission::GetFatalEnergyCheckLevels() const
00280 {
00281   // max energy non-conservation is mass of heavy nucleus
00282   return std::pair<G4double, G4double>(5*perCent,250*GeV);
00283 }

Generated on Mon May 27 17:48:47 2013 for Geant4 by  doxygen 1.4.7