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 #include "G4HadronElastic.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4ParticleTable.hh"
00036 #include "G4ParticleDefinition.hh"
00037 #include "G4IonTable.hh"
00038 #include "Randomize.hh"
00039 #include "G4Proton.hh"
00040 #include "G4Neutron.hh"
00041 #include "G4Deuteron.hh"
00042 #include "G4Alpha.hh"
00043 #include "G4Pow.hh"
00044
00045 G4HadronElastic::G4HadronElastic(const G4String& name)
00046 : G4HadronicInteraction(name)
00047 {
00048 SetMinEnergy( 0.0*GeV );
00049 SetMaxEnergy( 100.*TeV );
00050 lowestEnergyLimit= 1.e-6*eV;
00051
00052 theProton = G4Proton::Proton();
00053 theNeutron = G4Neutron::Neutron();
00054 theDeuteron = G4Deuteron::Deuteron();
00055 theAlpha = G4Alpha::Alpha();
00056
00057 }
00058
00059
00060 G4HadronElastic::~G4HadronElastic()
00061 {}
00062
00063
00064 void G4HadronElastic::Description() const
00065 {
00066 char* dirName = getenv("G4PhysListDocDir");
00067 if (dirName) {
00068 std::ofstream outFile;
00069 G4String outFileName = GetModelName() + ".html";
00070 G4String pathName = G4String(dirName) + "/" + outFileName;
00071 outFile.open(pathName);
00072 outFile << "<html>\n";
00073 outFile << "<head>\n";
00074
00075 outFile << "<title>Description of G4HadronElastic Model</title>\n";
00076 outFile << "</head>\n";
00077 outFile << "<body>\n";
00078
00079 outFile << "G4HadronElastic is a hadron-nucleus elastic scattering\n"
00080 << "model which uses the Gheisha two-exponential momentum\n"
00081 << "transfer parameterization. The model is fully relativistic\n"
00082 << "as opposed to the original Gheisha model which was not.\n"
00083 << "This model may be used for all long-lived hadrons at all\n"
00084 << "incident energies.\n";
00085
00086 outFile << "</body>\n";
00087 outFile << "</html>\n";
00088 outFile.close();
00089 }
00090 }
00091
00092
00093 G4HadFinalState* G4HadronElastic::ApplyYourself(
00094 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
00095 {
00096 theParticleChange.Clear();
00097
00098 const G4HadProjectile* aParticle = &aTrack;
00099 G4double ekin = aParticle->GetKineticEnergy();
00100 if(ekin <= lowestEnergyLimit) {
00101 theParticleChange.SetEnergyChange(ekin);
00102 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00103 return &theParticleChange;
00104 }
00105
00106 G4int A = targetNucleus.GetA_asInt();
00107 G4int Z = targetNucleus.GetZ_asInt();
00108
00109 G4double plab = aParticle->GetTotalMomentum();
00110
00111
00112 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00113 G4double m1 = theParticle->GetPDGMass();
00114
00115 if (verboseLevel>1) {
00116 G4cout << "G4HadronElastic: "
00117 << aParticle->GetDefinition()->GetParticleName()
00118 << " Plab(GeV/c)= " << plab/GeV
00119 << " Ekin(MeV) = " << ekin/MeV
00120 << " scattered off Z= " << Z
00121 << " A= " << A
00122 << G4endl;
00123 }
00124
00125 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z);
00126 G4LorentzVector lv1 = aParticle->Get4Momentum();
00127 G4LorentzVector lv(0.0,0.0,0.0,mass2);
00128 lv += lv1;
00129
00130 G4ThreeVector bst = lv.boostVector();
00131 lv1.boost(-bst);
00132
00133 G4ThreeVector p1 = lv1.vect();
00134 G4double momentumCMS = p1.mag();
00135 G4double tmax = 4.0*momentumCMS*momentumCMS;
00136
00137
00138 G4double t = SampleInvariantT(theParticle, plab, Z, A);
00139 G4double phi = G4UniformRand()*CLHEP::twopi;
00140 G4double cost = 1. - 2.0*t/tmax;
00141 G4double sint;
00142
00143
00144 if(cost > 1.0 || cost < -1.0) {
00145
00146 G4cout << "G4HadronElastic WARNING (1 - cost)= " << 1 - cost
00147 << " after scattering of "
00148 << aParticle->GetDefinition()->GetParticleName()
00149 << " p(GeV/c)= " << plab/GeV
00150 << " on an ion Z= " << Z << " A= " << A
00151 << G4endl;
00152
00153 cost = 1.0;
00154 sint = 0.0;
00155
00156
00157 } else {
00158 sint = std::sqrt((1.0-cost)*(1.0+cost));
00159 }
00160 if (verboseLevel>1) {
00161 G4cout << " t= " << t << " tmax(GeV^2)= " << tmax/(GeV*GeV)
00162 << " Pcms(GeV)= " << momentumCMS/GeV << " cos(t)=" << cost
00163 << " sin(t)=" << sint << G4endl;
00164 }
00165 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
00166 v1 *= momentumCMS;
00167 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),
00168 std::sqrt(momentumCMS*momentumCMS + m1*m1));
00169
00170 nlv1.boost(bst);
00171
00172 G4double eFinal = nlv1.e() - m1;
00173 if (verboseLevel > 1) {
00174 G4cout <<" m= " << m1 << " Efin(MeV)= " << eFinal
00175 << " Proj: 4-mom " << lv1 << " Final: " << nlv1
00176 << G4endl;
00177 }
00178 if(eFinal <= lowestEnergyLimit) {
00179 if(eFinal < 0.0 && verboseLevel > 0) {
00180 G4cout << "G4HadronElastic WARNING Efinal= " << eFinal
00181 << " after scattering of "
00182 << aParticle->GetDefinition()->GetParticleName()
00183 << " p(GeV/c)= " << plab/GeV
00184 << " on an ion Z= " << Z << " A= " << A
00185 << G4endl;
00186 }
00187 theParticleChange.SetEnergyChange(0.0);
00188 nlv1 = G4LorentzVector(0.0,0.0,0.0,m1);
00189
00190 } else {
00191 theParticleChange.SetMomentumChange(nlv1.vect().unit());
00192 theParticleChange.SetEnergyChange(eFinal);
00193 }
00194
00195 lv -= nlv1;
00196 G4double erec = lv.e() - mass2;
00197 if (verboseLevel > 1) {
00198 G4cout << "Recoil: " <<" m= " << mass2 << " Erec(MeV)= " << erec
00199 << " 4-mom: " << lv
00200 << G4endl;
00201 }
00202
00203 if(erec > GetRecoilEnergyThreshold()) {
00204 G4ParticleDefinition * theDef = 0;
00205 if(Z == 1 && A == 1) { theDef = theProton; }
00206 else if (Z == 1 && A == 2) { theDef = theDeuteron; }
00207 else if (Z == 1 && A == 3) { theDef = G4Triton::Triton(); }
00208 else if (Z == 2 && A == 3) { theDef = G4He3::He3(); }
00209 else if (Z == 2 && A == 4) { theDef = theAlpha; }
00210 else {
00211 theDef =
00212 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
00213 }
00214 G4DynamicParticle * aSec = new G4DynamicParticle(theDef, lv);
00215 theParticleChange.AddSecondary(aSec);
00216 } else if(erec > 0.0) {
00217 theParticleChange.SetLocalEnergyDeposit(erec);
00218 }
00219
00220 return &theParticleChange;
00221 }
00222
00223
00224 G4double
00225 G4HadronElastic::SampleInvariantT(const G4ParticleDefinition* p,
00226 G4double plab,
00227 G4int Z, G4int A)
00228 {
00229 static const G4double GeV2 = GeV*GeV;
00230 G4double momentumCMS = ComputeMomentumCMS(p,plab,Z,A);
00231 G4double tmax = 4.0*momentumCMS*momentumCMS/GeV2;
00232 G4double aa, bb, cc;
00233 G4double dd = 10.;
00234 G4Pow* g4pow = G4Pow::GetInstance();
00235 if (A <= 62) {
00236 bb = 14.5*g4pow->Z23(A);
00237 aa = g4pow->powZ(A, 1.63)/bb;
00238 cc = 1.4*g4pow->Z13(A)/dd;
00239 } else {
00240 bb = 60.*g4pow->Z13(A);
00241 aa = g4pow->powZ(A, 1.33)/bb;
00242 cc = 0.4*g4pow->powZ(A, 0.4)/dd;
00243 }
00244 G4double q1 = 1.0 - std::exp(-bb*tmax);
00245 G4double q2 = 1.0 - std::exp(-dd*tmax);
00246 G4double s1 = q1*aa;
00247 G4double s2 = q2*cc;
00248 if((s1 + s2)*G4UniformRand() < s2) {
00249 q1 = q2;
00250 bb = dd;
00251 }
00252 return -GeV2*std::log(1.0 - G4UniformRand()*q1)/bb;
00253 }
00254
