#include <G4LElastic.hh>
Inheritance diagram for G4LElastic:
Public Member Functions | |
| G4LElastic (const G4String &name="G4LElastic") | |
| ~G4LElastic () | |
| G4HadFinalState * | ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus) |
| virtual void | ModelDescription (std::ostream &outFile) const |
| virtual const std::pair< G4double, G4double > | GetFatalEnergyCheckLevels () const |
Definition at line 61 of file G4LElastic.hh.
| G4LElastic::G4LElastic | ( | const G4String & | name = "G4LElastic" |
) |
Definition at line 47 of file G4LElastic.cc.
References DBL_MAX, G4HadronicInteraction::SetMaxEnergy(), and G4HadronicInteraction::SetMinEnergy().
00048 :G4HadronicInteraction(name) 00049 { 00050 SetMinEnergy(0.0); 00051 SetMaxEnergy(DBL_MAX); 00052 }
| G4LElastic::~G4LElastic | ( | ) | [inline] |
| G4HadFinalState * G4LElastic::ApplyYourself | ( | const G4HadProjectile & | aTrack, | |
| G4Nucleus & | targetNucleus | |||
| ) | [virtual] |
Implements G4HadronicInteraction.
Definition at line 68 of file G4LElastic.cc.
References G4HadFinalState::AddSecondary(), G4AntiLambda::AntiLambda(), G4LightMedia::AntiLambdaExchange(), G4AntiNeutron::AntiNeutron(), G4LightMedia::AntiNeutronExchange(), G4AntiOmegaMinus::AntiOmegaMinus(), G4LightMedia::AntiOmegaMinusExchange(), G4AntiProton::AntiProton(), G4LightMedia::AntiProtonExchange(), G4AntiSigmaMinus::AntiSigmaMinus(), G4LightMedia::AntiSigmaMinusExchange(), G4AntiSigmaPlus::AntiSigmaPlus(), G4LightMedia::AntiSigmaPlusExchange(), G4AntiXiMinus::AntiXiMinus(), G4LightMedia::AntiXiMinusExchange(), G4AntiXiZero::AntiXiZero(), G4LightMedia::AntiXiZeroExchange(), G4HadFinalState::Clear(), G4ParticleTable::FindIon(), G4cout, G4endl, G4UniformRand, G4HadProjectile::Get4Momentum(), G4Nucleus::GetA_asInt(), G4HadProjectile::GetDefinition(), G4HadProjectile::GetKineticEnergy(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4HadProjectile::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4KaonMinus::KaonMinus(), G4LightMedia::KaonMinusExchange(), G4KaonPlus::KaonPlus(), G4LightMedia::KaonPlusExchange(), G4KaonZeroLong::KaonZeroLong(), G4LightMedia::KaonZeroLongExchange(), G4KaonZeroShort::KaonZeroShort(), G4LightMedia::KaonZeroShortExchange(), G4Lambda::Lambda(), G4LightMedia::LambdaExchange(), G4Neutron::Neutron(), G4LightMedia::NeutronExchange(), G4OmegaMinus::OmegaMinus(), G4LightMedia::OmegaMinusExchange(), G4PionMinus::PionMinus(), G4LightMedia::PionMinusExchange(), G4PionPlus::PionPlus(), G4LightMedia::PionPlusExchange(), G4Proton::Proton(), G4LightMedia::ProtonExchange(), G4HadFinalState::SetEnergyChange(), G4DynamicParticle::SetMomentum(), G4HadFinalState::SetMomentumChange(), G4DynamicParticle::SetMomentumDirection(), G4HadFinalState::SetStatusChange(), G4SigmaMinus::SigmaMinus(), G4LightMedia::SigmaMinusExchange(), G4SigmaPlus::SigmaPlus(), G4LightMedia::SigmaPlusExchange(), stopAndKill, G4HadronicInteraction::theParticleChange, G4HadronicInteraction::verboseLevel, G4XiMinus::XiMinus(), G4LightMedia::XiMinusExchange(), G4XiZero::XiZero(), and G4LightMedia::XiZeroExchange().
Referenced by G4NeutronHPorLElasticModel::ApplyYourself().
00069 { 00070 if(getenv("debug_LElastic")) verboseLevel = 5; 00071 theParticleChange.Clear(); 00072 const G4HadProjectile* aParticle = &aTrack; 00073 G4double atno2 = targetNucleus.GetA_asInt(); 00074 G4double zTarget = targetNucleus.GetZ_asInt(); 00075 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy()); 00076 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit()); 00077 00078 // Elastic scattering off Hydrogen 00079 00080 G4DynamicParticle* aSecondary = 0; 00081 if (atno2 < 1.5) { 00082 const G4ParticleDefinition* aParticleType = aParticle->GetDefinition(); 00083 if (aParticleType == G4PionPlus::PionPlus()) 00084 aSecondary = LightMedia.PionPlusExchange(aParticle, targetNucleus); 00085 else if (aParticleType == G4PionMinus::PionMinus()) 00086 aSecondary = LightMedia.PionMinusExchange(aParticle, targetNucleus); 00087 else if (aParticleType == G4KaonPlus::KaonPlus()) 00088 aSecondary = LightMedia.KaonPlusExchange(aParticle, targetNucleus); 00089 else if (aParticleType == G4KaonZeroShort::KaonZeroShort()) 00090 aSecondary = LightMedia.KaonZeroShortExchange(aParticle,targetNucleus); 00091 else if (aParticleType == G4KaonZeroLong::KaonZeroLong()) 00092 aSecondary = LightMedia.KaonZeroLongExchange(aParticle, targetNucleus); 00093 else if (aParticleType == G4KaonMinus::KaonMinus()) 00094 aSecondary = LightMedia.KaonMinusExchange(aParticle, targetNucleus); 00095 else if (aParticleType == G4Proton::Proton()) 00096 aSecondary = LightMedia.ProtonExchange(aParticle, targetNucleus); 00097 else if (aParticleType == G4AntiProton::AntiProton()) 00098 aSecondary = LightMedia.AntiProtonExchange(aParticle, targetNucleus); 00099 else if (aParticleType == G4Neutron::Neutron()) 00100 aSecondary = LightMedia.NeutronExchange(aParticle, targetNucleus); 00101 else if (aParticleType == G4AntiNeutron::AntiNeutron()) 00102 aSecondary = LightMedia.AntiNeutronExchange(aParticle, targetNucleus); 00103 else if (aParticleType == G4Lambda::Lambda()) 00104 aSecondary = LightMedia.LambdaExchange(aParticle, targetNucleus); 00105 else if (aParticleType == G4AntiLambda::AntiLambda()) 00106 aSecondary = LightMedia.AntiLambdaExchange(aParticle, targetNucleus); 00107 else if (aParticleType == G4SigmaPlus::SigmaPlus()) 00108 aSecondary = LightMedia.SigmaPlusExchange(aParticle, targetNucleus); 00109 else if (aParticleType == G4SigmaMinus::SigmaMinus()) 00110 aSecondary = LightMedia.SigmaMinusExchange(aParticle, targetNucleus); 00111 else if (aParticleType == G4AntiSigmaPlus::AntiSigmaPlus()) 00112 aSecondary = LightMedia.AntiSigmaPlusExchange(aParticle,targetNucleus); 00113 else if (aParticleType == G4AntiSigmaMinus::AntiSigmaMinus()) 00114 aSecondary= LightMedia.AntiSigmaMinusExchange(aParticle,targetNucleus); 00115 else if (aParticleType == G4XiZero::XiZero()) 00116 aSecondary = LightMedia.XiZeroExchange(aParticle, targetNucleus); 00117 else if (aParticleType == G4XiMinus::XiMinus()) 00118 aSecondary = LightMedia.XiMinusExchange(aParticle, targetNucleus); 00119 else if (aParticleType == G4AntiXiZero::AntiXiZero()) 00120 aSecondary = LightMedia.AntiXiZeroExchange(aParticle, targetNucleus); 00121 else if (aParticleType == G4AntiXiMinus::AntiXiMinus()) 00122 aSecondary = LightMedia.AntiXiMinusExchange(aParticle, targetNucleus); 00123 else if (aParticleType == G4OmegaMinus::OmegaMinus()) 00124 aSecondary = LightMedia.OmegaMinusExchange(aParticle, targetNucleus); 00125 else if (aParticleType == G4AntiOmegaMinus::AntiOmegaMinus()) 00126 aSecondary= LightMedia.AntiOmegaMinusExchange(aParticle,targetNucleus); 00127 else if (aParticleType == G4KaonPlus::KaonPlus()) 00128 aSecondary = LightMedia.KaonPlusExchange(aParticle, targetNucleus); 00129 } 00130 00131 // Has a charge or strangeness exchange occurred? 00132 if (aSecondary) { 00133 aSecondary->SetMomentum(aParticle->Get4Momentum().vect()); 00134 theParticleChange.SetStatusChange(stopAndKill); 00135 theParticleChange.AddSecondary(aSecondary); 00136 } 00137 // G4cout << "Entering elastic scattering 1"<<G4endl; 00138 00139 G4double p = aParticle->GetTotalMomentum()/GeV; 00140 if (verboseLevel > 1) 00141 G4cout << "G4LElastic::DoIt: Incident particle p=" << p << " GeV" << G4endl; 00142 00143 if (p < 0.01) return &theParticleChange; 00144 00145 // G4cout << "Entering elastic scattering 2"<<G4endl; 00146 // Compute the direction of elastic scattering. 00147 // It is planned to replace this code with a method based on 00148 // parameterized functions and a Monte Carlo method to invert the CDF. 00149 00150 G4double ran = G4UniformRand(); 00151 G4double aa, bb, cc, dd, rr; 00152 if (atno2 <= 62.) { 00153 aa = std::pow(atno2, 1.63); 00154 bb = 14.5*std::pow(atno2, 0.66); 00155 cc = 1.4*std::pow(atno2, 0.33); 00156 dd = 10.; 00157 } 00158 else { 00159 aa = std::pow(atno2, 1.33); 00160 bb = 60.*std::pow(atno2, 0.33); 00161 cc = 0.4*std::pow(atno2, 0.40); 00162 dd = 10.; 00163 } 00164 aa = aa/bb; 00165 cc = cc/dd; 00166 rr = (aa + cc)*ran; 00167 if (verboseLevel > 1) { 00168 G4cout << "DoIt: aa,bb,cc,dd,rr" << G4endl; 00169 G4cout << aa << " " << bb << " " << cc << " " << dd << " " << rr << G4endl; 00170 } 00171 G4double t1 = -std::log(ran)/bb; 00172 G4double t2 = -std::log(ran)/dd; 00173 if (verboseLevel > 1) { 00174 G4cout << "t1,Fctcos " << t1 << " " << Fctcos(t1, aa, bb, cc, dd, rr) << 00175 G4endl; 00176 G4cout << "t2,Fctcos " << t2 << " " << Fctcos(t2, aa, bb, cc, dd, rr) << 00177 G4endl; 00178 } 00179 G4double eps = 0.001; 00180 G4int ind1 = 10; 00181 G4double t; 00182 G4int ier1; 00183 ier1 = Rtmi(&t, t1, t2, eps, ind1, 00184 aa, bb, cc, dd, rr); 00185 if (verboseLevel > 1) { 00186 G4cout << "From Rtmi, ier1=" << ier1 << G4endl; 00187 G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << 00188 G4endl; 00189 } 00190 if (ier1 != 0) t = 0.25*(3.*t1 + t2); 00191 if (verboseLevel > 1) { 00192 G4cout << "t, Fctcos " << t << " " << Fctcos(t, aa, bb, cc, dd, rr) << 00193 G4endl; 00194 } 00195 G4double phi = G4UniformRand()*twopi; 00196 rr = 0.5*t/(p*p); 00197 if (rr > 1.) rr = 0.; 00198 if (verboseLevel > 1) 00199 G4cout << "rr=" << rr << G4endl; 00200 G4double cost = 1. - rr; 00201 G4double sint = std::sqrt(std::max(rr*(2. - rr), 0.)); 00202 if (sint == 0.) return &theParticleChange; 00203 // G4cout << "Entering elastic scattering 3"<<G4endl; 00204 if (verboseLevel > 1) 00205 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl; 00206 00207 // Scattered particle referred to axis of incident particle 00208 G4double mass1 = aParticle->GetDefinition()->GetPDGMass(); 00209 G4int Z=static_cast<G4int>(zTarget+.5); 00210 G4int A=static_cast<G4int>(atno2); 00211 if(G4UniformRand()<atno2-A) A++; 00212 //G4cout << " ion info "<<atno2 << " "<<A<<" "<<Z<<" "<<zTarget<<G4endl; 00213 G4double mass2 = 00214 G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z)->GetPDGMass(); 00215 00216 // non relativistic approximation 00217 G4double a = 1 + mass2/mass1; 00218 G4double b = -2.*p*cost; 00219 G4double c = p*p*(1 - mass2/mass1); 00220 G4double p1 = (-b+std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a); 00221 G4double px = p1*sint*std::sin(phi); 00222 G4double py = p1*sint*std::cos(phi); 00223 G4double pz = p1*cost; 00224 00225 // relativistic calculation 00226 G4double etot = std::sqrt(mass1*mass1+p*p)+mass2; 00227 a = etot*etot-p*p*cost*cost; 00228 b = 2*p*p*(mass1*cost*cost-etot); 00229 c = p*p*p*p*sint*sint; 00230 00231 G4double de = (-b-std::sqrt(std::max(0.0,b*b-4.*a*c)))/(2.*a); 00232 G4double e1 = std::sqrt(p*p+mass1*mass1)-de; 00233 G4double p12=e1*e1-mass1*mass1; 00234 p1 = std::sqrt(std::max(1.*eV*eV,p12)); 00235 px = p1*sint*std::sin(phi); 00236 py = p1*sint*std::cos(phi); 00237 pz = p1*cost; 00238 00239 if (verboseLevel > 1) 00240 { 00241 G4cout << "Relevant test "<<p<<" "<<p1<<" "<<cost<<" "<<de<<G4endl; 00242 G4cout << "p1/p = "<<p1/p<<" "<<mass1<<" "<<mass2<<" "<<a<<" "<<b<<" "<<c<<G4endl; 00243 G4cout << "rest = "<< b*b<<" "<<4.*a*c<<" "<<G4endl; 00244 G4cout << "make p1 = "<< p12<<" "<<e1*e1<<" "<<mass1*mass1<<" "<<G4endl; 00245 } 00246 // Incident particle 00247 G4double pxinc = p*aParticle->Get4Momentum().vect().unit().x(); 00248 G4double pyinc = p*aParticle->Get4Momentum().vect().unit().y(); 00249 G4double pzinc = p*aParticle->Get4Momentum().vect().unit().z(); 00250 if (verboseLevel > 1) { 00251 G4cout << "NOM SCAT " << px << " " << py << " " << pz << G4endl; 00252 G4cout << "INCIDENT " << pxinc << " " << pyinc << " " << pzinc << G4endl; 00253 } 00254 00255 // Transform scattered particle to reflect direction of incident particle 00256 G4double pxnew, pynew, pznew; 00257 Defs1(p, px, py, pz, pxinc, pyinc, pzinc, &pxnew, &pynew, &pznew); 00258 // Normalize: 00259 G4double pxre=pxinc-pxnew; 00260 G4double pyre=pyinc-pynew; 00261 G4double pzre=pzinc-pznew; 00262 G4ThreeVector it0(pxnew*GeV, pynew*GeV, pznew*GeV); 00263 if(p1>0) 00264 { 00265 pxnew = pxnew/p1; 00266 pynew = pynew/p1; 00267 pznew = pznew/p1; 00268 } 00269 else 00270 { 00271 //G4double pphi = 2*pi*G4UniformRand(); 00272 //G4double ccth = 2*G4UniformRand()-1; 00273 pxnew = 0;//std::sin(std::acos(ccth))*std::sin(pphi); 00274 pynew = 0;//std::sin(std::acos(ccth))*std::cos(phi); 00275 pznew = 1;//ccth; 00276 } 00277 if (verboseLevel > 1) { 00278 G4cout << "DoIt: returning new momentum vector" << G4endl; 00279 G4cout << "DoIt: "<<pxinc << " " << pyinc << " " << pzinc <<" "<<p<< G4endl; 00280 G4cout << "DoIt: "<<pxnew << " " << pynew << " " << pznew <<" "<<p<< G4endl; 00281 } 00282 00283 if (aSecondary) { 00284 aSecondary->SetMomentumDirection(pxnew, pynew, pznew); 00285 } else { 00286 try 00287 { 00288 theParticleChange.SetMomentumChange(pxnew, pynew, pznew); 00289 theParticleChange.SetEnergyChange(std::sqrt(mass1*mass1+it0.mag2())-mass1); 00290 } 00291 catch(G4HadronicException) 00292 { 00293 std::cerr << "GHADException originating from components of G4LElastic"<<std::cout; 00294 throw; 00295 } 00296 G4ParticleDefinition * theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z); 00297 G4ThreeVector it(pxre*GeV, pyre*GeV, pzre*GeV); 00298 G4DynamicParticle * aSec = 00299 new G4DynamicParticle(theDef, it.unit(), it.mag2()/(2.*mass2)); 00300 theParticleChange.AddSecondary(aSec); 00301 // G4cout << "Final check ###### "<<p<<" "<<it.mag()<<" "<<p1<<G4endl; 00302 } 00303 return &theParticleChange; 00304 }
Reimplemented from G4HadronicInteraction.
Definition at line 476 of file G4LElastic.cc.
00477 { 00478 return std::pair<G4double, G4double>(5*perCent,250*GeV); // max energy non-conservation is mass of heavy nucleus 00479 }
| void G4LElastic::ModelDescription | ( | std::ostream & | outFile | ) | const [virtual] |
Reimplemented from G4HadronicInteraction.
Definition at line 55 of file G4LElastic.cc.
00056 { 00057 outFile << "G4LElastic is one of the Low Energy Parameterized (LEP)\n" 00058 << "models used to implement elastic hadron scattering from nuclei.\n" 00059 << "It is a re-engineered version of the GHEISHA code of\n" 00060 << "H. Fesefeldt. It performs simplified two-body elastic\n" 00061 << "scattering for all long-lived hadronic projectiles by using\n" 00062 << "a two-exponential parameterization in momentum transfer.\n" 00063 << "It is valid for incident hadrons of all energies.\n"; 00064 }
1.4.7