#include <G4RPGTwoBody.hh>
Inheritance diagram for G4RPGTwoBody:
Public Member Functions | |
G4RPGTwoBody () | |
G4bool | ReactionStage (const G4HadProjectile *, G4ReactionProduct &, G4bool &, const G4DynamicParticle *, G4ReactionProduct &, G4bool &, const G4Nucleus &, G4ReactionProduct &, G4FastVector< G4ReactionProduct, 256 > &, G4int &, G4bool, G4ReactionProduct &) |
Definition at line 50 of file G4RPGTwoBody.hh.
G4RPGTwoBody::G4RPGTwoBody | ( | ) |
G4bool G4RPGTwoBody::ReactionStage | ( | const G4HadProjectile * | , | |
G4ReactionProduct & | , | |||
G4bool & | , | |||
const G4DynamicParticle * | , | |||
G4ReactionProduct & | , | |||
G4bool & | , | |||
const G4Nucleus & | , | |||
G4ReactionProduct & | , | |||
G4FastVector< G4ReactionProduct, 256 > & | , | |||
G4int & | , | |||
G4bool | , | |||
G4ReactionProduct & | ||||
) |
Reimplemented from G4RPGReaction.
Definition at line 44 of file G4RPGTwoBody.cc.
References G4RPGReaction::AddBlackTrackParticles(), G4RPGReaction::Defs1(), G4Poisson(), G4UniformRand, G4Nucleus::GetA_asInt(), G4ReactionProduct::GetDefinition(), G4Nucleus::GetDTABlackTrackEnergy(), G4RPGReaction::GetFinalStateNucleons(), G4ReactionProduct::GetKineticEnergy(), G4ReactionProduct::GetMass(), G4ReactionProduct::GetMomentum(), G4ParticleDefinition::GetParticleSubType(), G4ParticleDefinition::GetPDGMass(), G4Nucleus::GetPNBlackTrackEnergy(), G4ReactionProduct::GetTotalEnergy(), G4ReactionProduct::GetTotalMomentum(), G4Nucleus::GetZ_asInt(), G4ReactionProduct::Lorentz(), G4RPGReaction::normal(), G4InuclParticleNames::pp, G4ReactionProduct::SetKineticEnergy(), G4ReactionProduct::SetMass(), G4ReactionProduct::SetMomentum(), G4ReactionProduct::SetTOF(), and G4ReactionProduct::SetTotalEnergy().
Referenced by G4RPGInelastic::CalculateMomenta().
00056 { 00057 // 00058 // Derived from H. Fesefeldt's original FORTRAN code TWOB 00059 // 00060 // Generation of momenta for elastic and quasi-elastic 2 body reactions 00061 // 00062 // The simple formula ds/d|t| = s0* std::exp(-b*|t|) is used. 00063 // The b values are parametrizations from experimental data. 00064 // Unavailable values are taken from those of similar reactions. 00065 // 00066 00067 const G4double ekOriginal = modifiedOriginal.GetKineticEnergy()/GeV; 00068 const G4double etOriginal = modifiedOriginal.GetTotalEnergy()/GeV; 00069 const G4double mOriginal = modifiedOriginal.GetMass()/GeV; 00070 const G4double pOriginal = modifiedOriginal.GetMomentum().mag()/GeV; 00071 G4double currentMass = currentParticle.GetMass()/GeV; 00072 G4double targetMass = targetParticle.GetDefinition()->GetPDGMass()/GeV; 00073 00074 targetMass = targetParticle.GetMass()/GeV; 00075 const G4double atomicWeight = targetNucleus.GetA_asInt(); 00076 00077 G4double etCurrent = currentParticle.GetTotalEnergy()/GeV; 00078 G4double pCurrent = currentParticle.GetTotalMomentum()/GeV; 00079 00080 G4double cmEnergy = std::sqrt( currentMass*currentMass + 00081 targetMass*targetMass + 00082 2.0*targetMass*etCurrent ); // in GeV 00083 00084 if (cmEnergy < 0.01) { // 2-body scattering not possible 00085 targetParticle.SetMass( 0.0 ); // flag that the target particle doesn't exist 00086 00087 } else { 00088 // Projectile momentum in cm 00089 00090 G4double pf = targetMass*pCurrent/cmEnergy; 00091 00092 // 00093 // Set beam and target in centre of mass system 00094 // 00095 G4ReactionProduct pseudoParticle[3]; 00096 00097 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" || 00098 targetParticle.GetDefinition()->GetParticleSubType() == "pi") { 00099 00100 // G4double pf1 = pOriginal*mOriginal/std::sqrt(2.*mOriginal*(mOriginal+etOriginal)); 00101 00102 pseudoParticle[0].SetMass( targetMass*GeV ); 00103 pseudoParticle[0].SetTotalEnergy( etOriginal*GeV ); 00104 pseudoParticle[0].SetMomentum( 0.0, 0.0, pOriginal*GeV ); 00105 00106 pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 ); 00107 pseudoParticle[1].SetMass( mOriginal*GeV ); 00108 pseudoParticle[1].SetKineticEnergy( 0.0 ); 00109 00110 } else { 00111 pseudoParticle[0].SetMass( currentMass*GeV ); 00112 pseudoParticle[0].SetTotalEnergy( etCurrent*GeV ); 00113 pseudoParticle[0].SetMomentum( 0.0, 0.0, pCurrent*GeV ); 00114 00115 pseudoParticle[1].SetMomentum( 0.0, 0.0, 0.0 ); 00116 pseudoParticle[1].SetMass( targetMass*GeV ); 00117 pseudoParticle[1].SetKineticEnergy( 0.0 ); 00118 } 00119 // 00120 // Transform into center of mass system 00121 // 00122 pseudoParticle[2] = pseudoParticle[0] + pseudoParticle[1]; 00123 pseudoParticle[0].Lorentz( pseudoParticle[0], pseudoParticle[2] ); 00124 pseudoParticle[1].Lorentz( pseudoParticle[1], pseudoParticle[2] ); 00125 // 00126 // Set final state masses and energies in centre of mass system 00127 // 00128 currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV ); 00129 targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV ); 00130 00131 // 00132 // Calculate slope b for elastic scattering on proton/neutron 00133 // 00134 const G4double cb = 0.01; 00135 const G4double b1 = 4.225; 00136 const G4double b2 = 1.795; 00137 G4double b = std::max( cb, b1+b2*std::log(pOriginal) ); 00138 00139 // 00140 // Get cm scattering angle by sampling t from tmin to tmax 00141 // 00142 G4double btrang = b * 4.0 * pf * pseudoParticle[0].GetMomentum().mag()/GeV; 00143 G4double exindt = std::exp(-btrang) - 1.0; 00144 G4double costheta = 1.0 + 2*std::log( 1.0+G4UniformRand()*exindt ) /btrang; 00145 costheta = std::max(-1., std::min(1., costheta) ); 00146 G4double sintheta = std::sqrt((1.0-costheta)*(1.0+costheta)); 00147 G4double phi = twopi * G4UniformRand(); 00148 00149 // 00150 // Calculate final state momenta in centre of mass system 00151 // 00152 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" || 00153 targetParticle.GetDefinition()->GetParticleSubType() == "pi") { 00154 00155 currentParticle.SetMomentum( -pf*sintheta*std::cos(phi)*GeV, 00156 -pf*sintheta*std::sin(phi)*GeV, 00157 -pf*costheta*GeV ); 00158 } else { 00159 00160 currentParticle.SetMomentum( pf*sintheta*std::cos(phi)*GeV, 00161 pf*sintheta*std::sin(phi)*GeV, 00162 pf*costheta*GeV ); 00163 } 00164 targetParticle.SetMomentum( -currentParticle.GetMomentum() ); 00165 00166 // 00167 // Transform into lab system 00168 // 00169 currentParticle.Lorentz( currentParticle, pseudoParticle[1] ); 00170 targetParticle.Lorentz( targetParticle, pseudoParticle[1] ); 00171 00172 // Rotate final state particle vectors wrt incident momentum 00173 00174 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen ); 00175 00176 // Subtract binding energy 00177 00178 G4double pp, pp1, ekin; 00179 if( atomicWeight >= 1.5 ) 00180 { 00181 const G4double cfa = 0.025*((atomicWeight-1.)/120.)*std::exp(-(atomicWeight-1.)/120.); 00182 pp1 = currentParticle.GetMomentum().mag()/MeV; 00183 if( pp1 >= 1.0 ) 00184 { 00185 ekin = currentParticle.GetKineticEnergy()/MeV - cfa*(1.0+0.5*normal())*GeV; 00186 ekin = std::max( 0.0001*GeV, ekin ); 00187 currentParticle.SetKineticEnergy( ekin*MeV ); 00188 pp = currentParticle.GetTotalMomentum()/MeV; 00189 currentParticle.SetMomentum( currentParticle.GetMomentum() * (pp/pp1) ); 00190 } 00191 pp1 = targetParticle.GetMomentum().mag()/MeV; 00192 if( pp1 >= 1.0 ) 00193 { 00194 ekin = targetParticle.GetKineticEnergy()/MeV - cfa*(1.0+normal()/2.)*GeV; 00195 ekin = std::max( 0.0001*GeV, ekin ); 00196 targetParticle.SetKineticEnergy( ekin*MeV ); 00197 pp = targetParticle.GetTotalMomentum()/MeV; 00198 targetParticle.SetMomentum( targetParticle.GetMomentum() * (pp/pp1) ); 00199 } 00200 } 00201 } 00202 00203 // Get number of final state nucleons and nucleons remaining in 00204 // target nucleus 00205 00206 std::pair<G4int, G4int> finalStateNucleons = 00207 GetFinalStateNucleons(originalTarget, vec, vecLen); 00208 G4int protonsInFinalState = finalStateNucleons.first; 00209 G4int neutronsInFinalState = finalStateNucleons.second; 00210 00211 G4int PinNucleus = std::max(0, 00212 G4int(targetNucleus.GetZ_asInt()) - protonsInFinalState); 00213 G4int NinNucleus = std::max(0, 00214 G4int(targetNucleus.GetA_asInt()-targetNucleus.GetZ_asInt()) - neutronsInFinalState); 00215 00216 if( atomicWeight >= 1.5 ) 00217 { 00218 // Add black track particles 00219 // npnb: number of proton/neutron black track particles 00220 // ndta: number of deuterons, tritons, and alphas produced 00221 // epnb: kinetic energy available for proton/neutron black track 00222 // particles 00223 // edta: kinetic energy available for deuteron/triton/alpha particles 00224 00225 G4double epnb, edta; 00226 G4int npnb=0, ndta=0; 00227 00228 epnb = targetNucleus.GetPNBlackTrackEnergy(); // was enp1 in fortran code 00229 edta = targetNucleus.GetDTABlackTrackEnergy(); // was enp3 in fortran code 00230 const G4double pnCutOff = 0.0001; // GeV 00231 const G4double dtaCutOff = 0.0001; // GeV 00232 // const G4double kineticMinimum = 0.0001; 00233 // const G4double kineticFactor = -0.010; 00234 // G4double sprob = 0.0; // sprob = probability of self-absorption in heavy molecules 00235 if( epnb >= pnCutOff ) 00236 { 00237 npnb = G4Poisson( epnb/0.02 ); 00238 if( npnb > atomicWeight )npnb = G4int(atomicWeight); 00239 if( (epnb > pnCutOff) && (npnb <= 0) )npnb = 1; 00240 npnb = std::min( npnb, 127-vecLen ); 00241 } 00242 if( edta >= dtaCutOff ) 00243 { 00244 ndta = G4int(2.0 * std::log(atomicWeight)); 00245 ndta = std::min( ndta, 127-vecLen ); 00246 } 00247 00248 if (npnb == 0 && ndta == 0) npnb = 1; 00249 00250 AddBlackTrackParticles(epnb, npnb, edta, ndta, modifiedOriginal, 00251 PinNucleus, NinNucleus, targetNucleus, 00252 vec, vecLen); 00253 } 00254 00255 // 00256 // calculate time delay for nuclear reactions 00257 // 00258 if( (atomicWeight >= 1.5) && (atomicWeight <= 230.0) && (ekOriginal <= 0.2) ) 00259 currentParticle.SetTOF( 1.0-500.0*std::exp(-ekOriginal/0.04)*std::log(G4UniformRand()) ); 00260 else 00261 currentParticle.SetTOF( 1.0 ); 00262 00263 return true; 00264 }