G4LFission Class Reference

#include <G4LFission.hh>

Inheritance diagram for G4LFission:

G4HadronicInteraction

Public Member Functions

 G4LFission (const G4String &name="G4LFission")
 ~G4LFission ()
G4HadFinalStateApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
virtual void ModelDescription (std::ostream &outFile) const
virtual const std::pair< G4double,
G4double
GetFatalEnergyCheckLevels () const

Static Public Member Functions

static G4double Atomas (const G4double A, const G4double Z)

Detailed Description

Definition at line 65 of file G4LFission.hh.


Constructor & Destructor Documentation

G4LFission::G4LFission ( const G4String name = "G4LFission"  ) 

Definition at line 50 of file G4LFission.cc.

References DBL_MAX, G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().

00051  : G4HadronicInteraction(name)
00052 {
00053   init();
00054   SetMinEnergy(0.0*GeV);
00055   SetMaxEnergy(DBL_MAX);
00056 }

G4LFission::~G4LFission (  ) 

Definition at line 59 of file G4LFission.cc.

References G4HadFinalState::Clear(), and G4HadronicInteraction::theParticleChange.

00060 {
00061   theParticleChange.Clear();
00062 }


Member Function Documentation

G4HadFinalState * G4LFission::ApplyYourself ( const G4HadProjectile aTrack,
G4Nucleus targetNucleus 
) [virtual]

Implements G4HadronicInteraction.

Definition at line 94 of file G4LFission.cc.

References G4HadFinalState::AddSecondary(), Atomas(), G4HadFinalState::Clear(), G4cout, G4endl, G4UniformRand, G4Gamma::GammaDefinition(), G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4DynamicParticle::GetDefinition(), G4HadProjectile::GetDefinition(), G4HadProjectile::GetKineticEnergy(), G4HadSecondary::GetParticle(), G4ParticleDefinition::GetPDGCharge(), G4ParticleDefinition::GetPDGMass(), G4HadFinalState::GetSecondary(), G4DynamicParticle::GetTotalEnergy(), G4HadProjectile::GetTotalEnergy(), G4DynamicParticle::GetTotalMomentum(), G4HadProjectile::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4Neutron::NeutronDefinition(), G4InuclParticleNames::nn, G4InuclParticleNames::pp, G4DynamicParticle::SetKineticEnergy(), G4DynamicParticle::SetMomentumDirection(), G4HadFinalState::SetStatusChange(), stopAndKill, G4HadronicInteraction::theParticleChange, and G4HadronicInteraction::verboseLevel.

Referenced by G4NeutronHPorLFissionModel::ApplyYourself().

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 }

G4double G4LFission::Atomas ( const G4double  A,
const G4double  Z 
) [static]

Definition at line 243 of file G4LFission.cc.

References G4Alpha::AlphaDefinition(), G4Deuteron::DeuteronDefinition(), G4Electron::ElectronDefinition(), G4ParticleDefinition::GetPDGMass(), G4Neutron::NeutronDefinition(), and G4Proton::ProtonDefinition().

Referenced by ApplyYourself().

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 }

const std::pair< G4double, G4double > G4LFission::GetFatalEnergyCheckLevels (  )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 279 of file G4LFission.cc.

00280 {
00281   // max energy non-conservation is mass of heavy nucleus
00282   return std::pair<G4double, G4double>(5*perCent,250*GeV);
00283 }

void G4LFission::ModelDescription ( std::ostream &  outFile  )  const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 65 of file G4LFission.cc.

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 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:24 2013 for Geant4 by  doxygen 1.4.7