#include <G4HadronElastic.hh>
Inheritance diagram for G4HadronElastic:
Public Member Functions | |
G4HadronElastic (const G4String &name="hElasticLHEP") | |
virtual | ~G4HadronElastic () |
virtual G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
virtual G4double | SampleInvariantT (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) |
void | SetLowestEnergyLimit (G4double value) |
G4double | LowestEnergyLimit () const |
G4double | ComputeMomentumCMS (const G4ParticleDefinition *p, G4double plab, G4int Z, G4int A) |
virtual void | Description () const |
Definition at line 50 of file G4HadronElastic.hh.
G4HadronElastic::G4HadronElastic | ( | const G4String & | name = "hElasticLHEP" |
) |
Definition at line 45 of file G4HadronElastic.cc.
References G4Alpha::Alpha(), G4Deuteron::Deuteron(), G4Neutron::Neutron(), G4Proton::Proton(), G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().
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 //Description(); 00057 }
G4HadronElastic::~G4HadronElastic | ( | ) | [virtual] |
G4HadFinalState * G4HadronElastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
G4Nucleus & | targetNucleus | |||
) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 93 of file G4HadronElastic.cc.
References G4HadFinalState::AddSecondary(), G4HadFinalState::Clear(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4HadProjectile::GetKineticEnergy(), G4NucleiProperties::GetNuclearMass(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4HadronicInteraction::GetRecoilEnergyThreshold(), G4HadProjectile::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4He3::He3(), SampleInvariantT(), G4HadFinalState::SetEnergyChange(), G4HadFinalState::SetLocalEnergyDeposit(), G4HadFinalState::SetMomentumChange(), G4HadronicInteraction::theParticleChange, G4Triton::Triton(), and G4HadronicInteraction::verboseLevel.
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 // Scattered particle referred to axis of incident particle 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 // Sampling in CM system 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 // problem in sampling 00144 if(cost > 1.0 || cost < -1.0) { 00145 //if(verboseLevel > 0) { 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 // normal situation 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 }
G4double G4HadronElastic::ComputeMomentumCMS | ( | const G4ParticleDefinition * | p, | |
G4double | plab, | |||
G4int | Z, | |||
G4int | A | |||
) | [inline] |
Definition at line 98 of file G4HadronElastic.hh.
References G4NucleiProperties::GetNuclearMass(), and G4ParticleDefinition::GetPDGMass().
Referenced by SampleInvariantT().
00100 { 00101 G4double m1 = p->GetPDGMass(); 00102 G4double m12= m1*m1; 00103 G4double mass2 = G4NucleiProperties::GetNuclearMass(A, Z); 00104 return plab*mass2/std::sqrt(m12 + mass2*mass2 + 2.*mass2*std::sqrt(m12 + plab*plab)); 00105 }
void G4HadronElastic::Description | ( | ) | const [virtual] |
Reimplemented in G4CHIPSElastic, G4ChipsElasticModel, and G4ElasticHadrNucleusHE.
Definition at line 64 of file G4HadronElastic.cc.
References G4HadronicInteraction::GetModelName().
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 }
G4double G4HadronElastic::LowestEnergyLimit | ( | ) | const [inline] |
G4double G4HadronElastic::SampleInvariantT | ( | const G4ParticleDefinition * | p, | |
G4double | plab, | |||
G4int | Z, | |||
G4int | A | |||
) | [virtual] |
Reimplemented from G4HadronicInteraction.
Reimplemented in G4AntiNuclElastic, G4CHIPSElastic, G4ChipsElasticModel, G4DiffuseElastic, G4ElasticHadrNucleusHE, and G4NuclNuclDiffuseElastic.
Definition at line 225 of file G4HadronElastic.cc.
References ComputeMomentumCMS(), G4UniformRand, G4Pow::GetInstance(), G4Pow::powZ(), G4Pow::Z13(), and G4Pow::Z23().
Referenced by ApplyYourself(), G4ChipsElasticModel::SampleInvariantT(), and G4CHIPSElastic::SampleInvariantT().
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 }
void G4HadronElastic::SetLowestEnergyLimit | ( | G4double | value | ) | [inline] |
Reimplemented in G4DiffuseElastic, and G4NuclNuclDiffuseElastic.
Definition at line 87 of file G4HadronElastic.hh.