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 "G4LEDeuteronInelastic.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "Randomize.hh"
00034 #include "G4Electron.hh"
00035
00036 void G4LEDeuteronInelastic::ModelDescription(std::ostream& outFile) const
00037 {
00038 outFile << "G4LEDeuteronInelastic is one of the Low Energy Parameterized\n"
00039 << "(LEP) models used to implement inelastic deuteron 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 deuterons with initial energies between 0 and 10\n"
00046 << "TeV.\n";
00047 }
00048
00049
00050 G4HadFinalState*
00051 G4LEDeuteronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00052 G4Nucleus& targetNucleus)
00053 {
00054 theParticleChange.Clear();
00055 const G4HadProjectile* originalIncident = &aTrack;
00056
00057 if (verboseLevel > 1) {
00058 const G4Material *targetMaterial = aTrack.GetMaterial();
00059 G4cout << "G4LEDeuteronInelastic::ApplyYourself called" << G4endl;
00060 G4cout << "kinetic energy = " << originalIncident->GetKineticEnergy()/MeV << "MeV, ";
00061 G4cout << "target material = " << targetMaterial->GetName() << ", ";
00062 }
00063
00064
00065 if (originalIncident->GetKineticEnergy()/MeV > 100. ||
00066 originalIncident->GetKineticEnergy() <= 0.1*MeV) {
00067 theParticleChange.SetStatusChange(isAlive);
00068 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00069 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00070 return &theParticleChange;
00071 }
00072
00073 G4double A = targetNucleus.GetA_asInt();
00074 G4double Z = targetNucleus.GetZ_asInt();
00075 G4double theAtomicMass = targetNucleus.AtomicMass(A, Z);
00076 G4double massVec[9];
00077 massVec[0] = targetNucleus.AtomicMass( A+2.0, Z+1.0 );
00078 massVec[1] = targetNucleus.AtomicMass( A+1.0, Z+1.0 );
00079 massVec[2] = targetNucleus.AtomicMass( A+1.0, Z );
00080 massVec[3] = theAtomicMass;
00081 massVec[4] = 0.;
00082 if (A > 1.0 && A-1.0 > Z)
00083 massVec[4] = targetNucleus.AtomicMass(A-1.0, Z);
00084 massVec[5] = 0.;
00085 if (A > 2.0 && Z > 1.0 && A-2.0 > Z-1.0)
00086 massVec[5] = targetNucleus.AtomicMass(A-2.0, Z-1.0);
00087 massVec[6] = 0.;
00088 if (A > Z+1.0)
00089 massVec[6] = targetNucleus.AtomicMass(A, Z+1.0);
00090 massVec[7] = massVec[3];
00091 massVec[8] = 0.;
00092 if (Z > 1.0) massVec[8] = targetNucleus.AtomicMass(A,Z-1.0);
00093
00094 G4FastVector<G4ReactionProduct,4> vec;
00095 G4int vecLen = 0;
00096 vec.Initialize( 0 );
00097
00098 theReactionDynamics.NuclearReaction(vec, vecLen, originalIncident,
00099 targetNucleus, theAtomicMass, massVec);
00100
00101 G4double p = vec[0]->GetMomentum().mag();
00102 theParticleChange.SetMomentumChange( vec[0]->GetMomentum() * (1.0/p) );
00103 theParticleChange.SetEnergyChange( vec[0]->GetKineticEnergy() );
00104 delete vec[0];
00105
00106 if (vecLen <= 1)
00107 {
00108 theParticleChange.SetStatusChange(isAlive);
00109 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00110 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00111 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00112 return &theParticleChange;
00113 }
00114
00115 G4DynamicParticle* pd;
00116 for (G4int i=1; i<vecLen; ++i) {
00117 pd = new G4DynamicParticle();
00118 pd->SetDefinition( vec[i]->GetDefinition() );
00119 pd->SetMomentum( vec[i]->GetMomentum() );
00120 theParticleChange.AddSecondary( pd );
00121 delete vec[i];
00122 }
00123
00124 if (isotopeProduction) DoIsotopeCounting(originalIncident, targetNucleus);
00125 return &theParticleChange;
00126 }
00127
00128
00129