00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
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
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
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
00132
00133 G4double p;
00134 if (Z == 1 && N == 1) {
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
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
00198 return std::pair<G4double, G4double>(5*perCent,250*GeV);
00199 }