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 #include <iostream>
00030 #include <signal.h>
00031
00032 #include "G4RPGTwoBody.hh"
00033 #include "Randomize.hh"
00034 #include "G4PhysicalConstants.hh"
00035 #include "G4SystemOfUnits.hh"
00036 #include "G4Poisson.hh"
00037 #include "G4HadReentrentException.hh"
00038
00039 G4RPGTwoBody::G4RPGTwoBody()
00040 : G4RPGReaction() {}
00041
00042
00043 G4bool G4RPGTwoBody::
00044 ReactionStage(const G4HadProjectile* ,
00045 G4ReactionProduct& modifiedOriginal,
00046 G4bool& ,
00047 const G4DynamicParticle* originalTarget,
00048 G4ReactionProduct& targetParticle,
00049 G4bool& ,
00050 const G4Nucleus& targetNucleus,
00051 G4ReactionProduct& currentParticle,
00052 G4FastVector<G4ReactionProduct,256>& vec,
00053 G4int& vecLen,
00054 G4bool ,
00055 G4ReactionProduct& )
00056 {
00057
00058
00059
00060
00061
00062
00063
00064
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 );
00083
00084 if (cmEnergy < 0.01) {
00085 targetParticle.SetMass( 0.0 );
00086
00087 } else {
00088
00089
00090 G4double pf = targetMass*pCurrent/cmEnergy;
00091
00092
00093
00094
00095 G4ReactionProduct pseudoParticle[3];
00096
00097 if (targetParticle.GetDefinition()->GetParticleSubType() == "kaon" ||
00098 targetParticle.GetDefinition()->GetParticleSubType() == "pi") {
00099
00100
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
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
00127
00128 currentParticle.SetTotalEnergy( std::sqrt(pf*pf+currentMass*currentMass)*GeV );
00129 targetParticle.SetTotalEnergy( std::sqrt(pf*pf+targetMass*targetMass)*GeV );
00130
00131
00132
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
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
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
00168
00169 currentParticle.Lorentz( currentParticle, pseudoParticle[1] );
00170 targetParticle.Lorentz( targetParticle, pseudoParticle[1] );
00171
00172
00173
00174 Defs1( modifiedOriginal, currentParticle, targetParticle, vec, vecLen );
00175
00176
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
00204
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
00219
00220
00221
00222
00223
00224
00225 G4double epnb, edta;
00226 G4int npnb=0, ndta=0;
00227
00228 epnb = targetNucleus.GetPNBlackTrackEnergy();
00229 edta = targetNucleus.GetDTABlackTrackEnergy();
00230 const G4double pnCutOff = 0.0001;
00231 const G4double dtaCutOff = 0.0001;
00232
00233
00234
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
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 }
00265
00266