G4LCapture.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 //
00027 // $Id$
00028 //
00029 //
00030 // G4 Model: Low-energy Neutron Capture
00031 // F.W. Jones, TRIUMF, 03-DEC-96
00032 // 
00033 // This is a prototype of a low-energy neutron capture process.
00034 // Currently it is based on the GHEISHA routine CAPTUR,
00035 // and conforms fairly closely to the original Fortran.
00036 //
00037 // HPW Capture using models now. the code comes from the
00038 // original G4LCapture class.
00039 //
00040 // 25-JUN-98 FWJ: replaced missing Initialize for ParticleChange.
00041 //
00042 
00043 #include <iostream>
00044 
00045 #include "globals.hh"
00046 #include "Randomize.hh"
00047 #include "G4PhysicalConstants.hh"
00048 #include "G4SystemOfUnits.hh"
00049 #include "G4LCapture.hh"
00050 
00051 G4LCapture::G4LCapture(const G4String& name)
00052  : G4HadronicInteraction(name)
00053 {   
00054   SetMinEnergy(0.0*GeV);
00055   SetMaxEnergy(DBL_MAX);
00056   // Description();
00057 }
00058 
00059 
00060 G4LCapture::~G4LCapture()
00061 {
00062   theParticleChange.Clear();
00063 }
00064 
00065 
00066 void G4LCapture::ModelDescription(std::ostream& outFile) const
00067 {
00068   outFile << "G4LCapture is one of the Low Energy Parameterized\n"
00069           << "(LEP) models used to implement neutron capture on nuclei.\n"
00070           << "It is a re-engineered version of the GHEISHA code of\n"
00071           << "H. Fesefeldt which simply adds the neutron mass and energy\n"
00072           << "to the target nucleus, and emits gammas isotropically as\n"
00073           << "long as there is sufficient excitation energy in the\n"
00074           << "daughter nucleus.  The model is applicable to all incident\n"
00075           << "neutron energies.\n";
00076 }
00077 
00078 
00079 G4HadFinalState*
00080 G4LCapture::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
00081 {
00082   theParticleChange.Clear();
00083   theParticleChange.SetStatusChange(stopAndKill);
00084 
00085   G4double N = targetNucleus.GetA_asInt();
00086   G4double Z = targetNucleus.GetZ_asInt();
00087 
00088   const G4LorentzVector theMom = aTrack.Get4Momentum();
00089   G4double P = theMom.vect().mag()/GeV;
00090   G4double Px = theMom.vect().x()/GeV;
00091   G4double Py = theMom.vect().y()/GeV;
00092   G4double Pz = theMom.vect().z()/GeV;
00093   G4double E = theMom.e()/GeV;
00094   G4double E0 = aTrack.GetDefinition()->GetPDGMass()/GeV;
00095   G4double Q = aTrack.GetDefinition()->GetPDGCharge();
00096   if (verboseLevel > 1) {
00097     G4cout << "G4LCapture:ApplyYourself: incident particle:" << G4endl;
00098     G4cout << "P      " << P << " GeV/c" << G4endl;
00099     G4cout << "Px     " << Px << " GeV/c" << G4endl;
00100     G4cout << "Py     " << Py << " GeV/c" << G4endl;
00101     G4cout << "Pz     " << Pz << " GeV/c" << G4endl;
00102     G4cout << "E      " << E << " GeV" << G4endl;
00103     G4cout << "mass   " << E0 << " GeV" << G4endl;
00104     G4cout << "charge " << Q << G4endl;
00105   }
00106 
00107   // GHEISHA ADD operation to get total energy, mass, charge:
00108 
00109   if (verboseLevel > 1) {
00110     G4cout << "G4LCapture:ApplyYourself: material:" << G4endl;
00111     G4cout << "A      " << N << G4endl;
00112     G4cout << "Z   " << Z  << G4endl;
00113     G4cout << "atomic mass " << 
00114     Atomas(N, Z) << "GeV" << G4endl;
00115   }
00116   E = E + Atomas(N, Z);
00117   G4double E02 = E*E - P*P;
00118   E0 = std::sqrt(std::abs(E02));
00119   if (E02 < 0) E0 = -E0;
00120   Q = Q + Z;
00121   if (verboseLevel > 1) {
00122     G4cout << "G4LCapture:ApplyYourself: total:" << G4endl;
00123     G4cout << "E      " << E << " GeV" << G4endl;
00124     G4cout << "mass   " << E0 << " GeV" << G4endl;
00125     G4cout << "charge " << Q << G4endl;
00126   }
00127   Px = -Px;
00128   Py = -Py;
00129   Pz = -Pz;
00130 
00131   // Make a gamma...
00132 
00133   G4double p;
00134   if (Z == 1 && N == 1) { // special case for hydrogen
00135     p = 0.0022;
00136   } else {
00137     G4double ran = G4RandGauss::shoot();
00138     p = 0.0065 + ran*0.0010;
00139   }
00140 
00141   G4double ran1 = G4UniformRand();
00142   G4double ran2 = G4UniformRand();
00143   G4double cost = -1. + 2.*ran1;
00144   G4double sint = std::sqrt(std::abs(1. - cost*cost));
00145   G4double phi = ran2*twopi;
00146 
00147   G4double px = p*sint*std::sin(phi);
00148   G4double py = p*sint*std::cos(phi);
00149   G4double pz = p*cost;
00150   G4double e = p;
00151 
00152   G4double a = px*Px + py*Py + pz*Pz;
00153   a = (a/(E + E0) - e)/E0;
00154 
00155   px = px + a*Px;
00156   py = py + a*Py;
00157   pz = pz + a*Pz;
00158 
00159   G4DynamicParticle* aGamma;
00160   aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(),
00161                                  G4ThreeVector(px*GeV, py*GeV, pz*GeV));
00162   theParticleChange.AddSecondary(aGamma);
00163 
00164   // Make another gamma if there is sufficient energy left over...
00165 
00166   G4double xp = 0.008 - p;
00167   if (xp > 0.) {
00168     if (Z > 1 || N > 1) { 
00169       ran1 = G4UniformRand();
00170       ran2 = G4UniformRand();
00171       cost = -1. + 2.*ran1;
00172       sint = std::sqrt(std::abs(1. - cost*cost));
00173       phi = ran2*twopi;
00174 
00175       px = xp*sint*std::sin(phi);
00176       py = xp*sint*std::cos(phi);
00177       pz = xp*cost;
00178       e = xp;
00179 
00180       a = px*Px + py*Py + pz*Pz;
00181       a = (a/(E + E0) - e)/E0;
00182 
00183       px = px + a*Px;
00184       py = py + a*Py;
00185       pz = pz + a*Pz;
00186 
00187       aGamma = new G4DynamicParticle(G4Gamma::GammaDefinition(),
00188                                      G4ThreeVector(px*GeV, py*GeV, pz*GeV));
00189       theParticleChange.AddSecondary(aGamma);
00190     }
00191   }
00192   return &theParticleChange;
00193 }
00194 
00195 const std::pair<G4double, G4double> G4LCapture::GetFatalEnergyCheckLevels() const
00196 {
00197         // max energy non-conservation is mass of heavy nucleus
00198         return std::pair<G4double, G4double>(5*perCent,250*GeV);
00199 }

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