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 #include "G4LETritonInelastic.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "Randomize.hh"
00034 #include "G4Electron.hh"
00035
00036 void G4LETritonInelastic::ModelDescription(std::ostream& outFile) const
00037 {
00038 outFile << "G4LETritonInelastic is one of the Low Energy Parameterized\n"
00039 << "(LEP) models used to implement inelastic triton scattering\n"
00040 << "from nuclei. It is a re-engineered version of the GHEISHA\n"
00041 << "code of H. Fesefeldt. It divides the initial collision\n"
00042 << "products into backward- and forward-going clusters which are\n"
00043 << "then decayed into final state hadrons. The model does not\n"
00044 << "conserve energy on an event-by-event basis. It may be\n"
00045 << "applied to tritons with initial energies between 0 and 25\n"
00046 << "GeV.\n";
00047 }
00048
00049
00050 G4HadFinalState*
00051 G4LETritonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00052 G4Nucleus& targetNucleus)
00053 {
00054 G4bool triton_debug = false;
00055 if (getenv("TritonLEDebug")) triton_debug = true;
00056 theParticleChange.Clear();
00057 const G4HadProjectile *originalIncident = &aTrack;
00058 if (triton_debug) G4cout << "entering LETritonInelastic "
00059 << originalIncident->GetKineticEnergy() << G4endl;
00060 if (originalIncident->GetKineticEnergy() <= 0.1*MeV) {
00061 theParticleChange.SetStatusChange(isAlive);
00062 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00063 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00064 return &theParticleChange;
00065 }
00066
00067 if (verboseLevel > 1) {
00068 const G4Material *targetMaterial = aTrack.GetMaterial();
00069 G4cout << "G4LETritonInelastic::ApplyYourself called" << G4endl;
00070 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00071 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00072 }
00073
00074 if (triton_debug) G4cout << "running LETritonInelastic 1" << G4endl;
00075
00076
00077 if (originalIncident->GetKineticEnergy()/MeV > 100. ||
00078 originalIncident->GetKineticEnergy() <= 0.) {
00079 theParticleChange.SetStatusChange(isAlive);
00080 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00081 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00082 return &theParticleChange;
00083 }
00084
00085 if (triton_debug) G4cout << "running LETritonInelastic 2" << G4endl;
00086
00087 G4double A = targetNucleus.GetA_asInt();
00088 G4double Z = targetNucleus.GetZ_asInt();
00089 G4double theAtomicMass = targetNucleus.AtomicMass(A, Z);
00090 G4double massVec[9];
00091 massVec[0] = targetNucleus.AtomicMass( A+3.0, Z+1.0 );
00092 massVec[1] = targetNucleus.AtomicMass( A+2.0, Z+1.0 );
00093 massVec[2] = targetNucleus.AtomicMass( A+2.0, Z );
00094 massVec[3] = targetNucleus.AtomicMass( A+1.0, Z );
00095 massVec[4] = theAtomicMass;
00096 massVec[5] = massVec[3];
00097 if (A > 1.0 && Z > 1.0) massVec[5] = targetNucleus.AtomicMass(A-1.0, Z-1.0);
00098 massVec[6] = targetNucleus.AtomicMass(A+1.0, Z+1.0);
00099 massVec[7] = massVec[3];
00100 massVec[8] = massVec[2];
00101 if (Z > 1.0) massVec[8] = targetNucleus.AtomicMass(A+1.0, Z-1.0);
00102
00103 G4FastVector<G4ReactionProduct,4> vec;
00104 G4int vecLen = 0;
00105 vec.Initialize(0);
00106
00107 if (triton_debug) G4cout << "running LETritonInelastic 3" << G4endl;
00108 theReactionDynamics.NuclearReaction(vec, vecLen, originalIncident,
00109 targetNucleus, theAtomicMass, massVec);
00110 if (triton_debug) G4cout << "running LETritonInelastic 4" << G4endl;
00111
00112 G4double p = vec[0]->GetMomentum().mag();
00113 theParticleChange.SetMomentumChange( vec[0]->GetMomentum()*(1./p) );
00114 theParticleChange.SetEnergyChange( vec[0]->GetKineticEnergy() );
00115 delete vec[0];
00116
00117 if (vecLen <= 1) {
00118 theParticleChange.SetStatusChange(isAlive);
00119 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00120 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00121 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00122 return &theParticleChange;
00123 }
00124
00125 G4DynamicParticle *pd;
00126 for (G4int i = 1; i < vecLen; ++i) {
00127 pd = new G4DynamicParticle();
00128 pd->SetDefinition( vec[i]->GetDefinition() );
00129 pd->SetMomentum( vec[i]->GetMomentum() );
00130 theParticleChange.AddSecondary( pd );
00131 delete vec[i];
00132 }
00133
00134 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00135
00136 if (triton_debug) G4cout << "leaving LETritonInelastic" << G4endl;
00137 return &theParticleChange;
00138 }
00139
00140
00141