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 "G4RPGNeutronInelastic.hh"
00030 #include "G4PhysicalConstants.hh"
00031 #include "G4SystemOfUnits.hh"
00032 #include "Randomize.hh"
00033
00034 G4HadFinalState*
00035 G4RPGNeutronInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00036 G4Nucleus& targetNucleus)
00037 {
00038 theParticleChange.Clear();
00039 const G4HadProjectile* originalIncident = &aTrack;
00040
00041
00042 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00043
00044 G4ReactionProduct modifiedOriginal;
00045 modifiedOriginal = *originalIncident;
00046 G4ReactionProduct targetParticle;
00047 targetParticle = *originalTarget;
00048 if( originalIncident->GetKineticEnergy()/GeV < 0.01 + 2.*G4UniformRand()/9. )
00049 {
00050 SlowNeutron(originalIncident,modifiedOriginal,targetParticle,targetNucleus );
00051 delete originalTarget;
00052 return &theParticleChange;
00053 }
00054
00055
00056
00057 G4double ek = originalIncident->GetKineticEnergy()/MeV;
00058 G4double amas = originalIncident->GetDefinition()->GetPDGMass()/MeV;
00059
00060 G4double tkin = targetNucleus.Cinema( ek );
00061 ek += tkin;
00062 modifiedOriginal.SetKineticEnergy( ek*MeV );
00063 G4double et = ek + amas;
00064 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00065 G4double pp = modifiedOriginal.GetMomentum().mag()/MeV;
00066 if( pp > 0.0 )
00067 {
00068 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00069 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00070 }
00071
00072
00073
00074 tkin = targetNucleus.EvaporationEffects( ek );
00075 ek -= tkin;
00076 modifiedOriginal.SetKineticEnergy(ek);
00077 et = ek + amas;
00078 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00079 pp = modifiedOriginal.GetMomentum().mag();
00080 if( pp > 0.0 )
00081 {
00082 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00083 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00084 }
00085 const G4double cutOff = 0.1;
00086 if( modifiedOriginal.GetKineticEnergy()/MeV <= cutOff )
00087 {
00088 SlowNeutron( originalIncident, modifiedOriginal, targetParticle, targetNucleus );
00089 delete originalTarget;
00090 return &theParticleChange;
00091 }
00092
00093 G4ReactionProduct currentParticle = modifiedOriginal;
00094 currentParticle.SetSide( 1 );
00095 targetParticle.SetSide( -1 );
00096 G4bool incidentHasChanged = false;
00097 G4bool targetHasChanged = false;
00098 G4bool quasiElastic = false;
00099 G4FastVector<G4ReactionProduct,256> vec;
00100 G4int vecLen = 0;
00101 vec.Initialize( 0 );
00102
00103 InitialCollision(vec, vecLen, currentParticle, targetParticle,
00104 incidentHasChanged, targetHasChanged);
00105
00106 CalculateMomenta(vec, vecLen,
00107 originalIncident, originalTarget, modifiedOriginal,
00108 targetNucleus, currentParticle, targetParticle,
00109 incidentHasChanged, targetHasChanged, quasiElastic);
00110
00111 SetUpChange(vec, vecLen,
00112 currentParticle, targetParticle,
00113 incidentHasChanged);
00114
00115 delete originalTarget;
00116 return &theParticleChange;
00117 }
00118
00119 void
00120 G4RPGNeutronInelastic::SlowNeutron(const G4HadProjectile* originalIncident,
00121 G4ReactionProduct& modifiedOriginal,
00122 G4ReactionProduct& targetParticle,
00123 G4Nucleus& targetNucleus)
00124 {
00125 const G4double A = targetNucleus.GetA_asInt();
00126 const G4double Z = targetNucleus.GetZ_asInt();
00127
00128 G4double currentKinetic = modifiedOriginal.GetKineticEnergy()/MeV;
00129 G4double currentMass = modifiedOriginal.GetMass()/MeV;
00130 if( A < 1.5 )
00131 {
00132
00133
00134
00135
00136
00137 G4double cost1, eka = 0.0;
00138 while (eka <= 0.0)
00139 {
00140 cost1 = -1.0 + 2.0*G4UniformRand();
00141 eka = 1.0 + 2.0*cost1*A + A*A;
00142 }
00143 G4double cost = std::min( 1.0, std::max( -1.0, (A*cost1+1.0)/std::sqrt(eka) ) );
00144 eka /= (1.0+A)*(1.0+A);
00145 G4double ek = currentKinetic*MeV/GeV;
00146 G4double amas = currentMass*MeV/GeV;
00147 ek *= eka;
00148 G4double en = ek + amas;
00149 G4double p = std::sqrt(std::abs(en*en-amas*amas));
00150 G4double sint = std::sqrt(std::abs(1.0-cost*cost));
00151 G4double phi = G4UniformRand()*twopi;
00152 G4double px = sint*std::sin(phi);
00153 G4double py = sint*std::cos(phi);
00154 G4double pz = cost;
00155 targetParticle.SetMomentum( px*GeV, py*GeV, pz*GeV );
00156 G4double pxO = originalIncident->Get4Momentum().x()/GeV;
00157 G4double pyO = originalIncident->Get4Momentum().y()/GeV;
00158 G4double pzO = originalIncident->Get4Momentum().z()/GeV;
00159 G4double ptO = pxO*pxO + pyO+pyO;
00160 if( ptO > 0.0 )
00161 {
00162 G4double pO = std::sqrt(pxO*pxO+pyO*pyO+pzO*pzO);
00163 cost = pzO/pO;
00164 sint = 0.5*(std::sqrt(std::abs((1.0-cost)*(1.0+cost)))+std::sqrt(ptO)/pO);
00165 G4double ph = pi/2.0;
00166 if( pyO < 0.0 )ph = ph*1.5;
00167 if( std::abs(pxO) > 0.000001 )ph = std::atan2(pyO,pxO);
00168 G4double cosp = std::cos(ph);
00169 G4double sinp = std::sin(ph);
00170 px = cost*cosp*px - sinp*py+sint*cosp*pz;
00171 py = cost*sinp*px + cosp*py+sint*sinp*pz;
00172 pz = -sint*px + cost*pz;
00173 }
00174 else
00175 {
00176 if( pz < 0.0 )pz *= -1.0;
00177 }
00178 G4double pu = std::sqrt(px*px+py*py+pz*pz);
00179 modifiedOriginal.SetMomentum( targetParticle.GetMomentum() * (p/pu) );
00180 modifiedOriginal.SetKineticEnergy( ek*GeV );
00181
00182 targetParticle.SetMomentum(
00183 originalIncident->Get4Momentum().vect() - modifiedOriginal.GetMomentum() );
00184 G4double pp = targetParticle.GetMomentum().mag();
00185 G4double tarmas = targetParticle.GetMass();
00186 targetParticle.SetTotalEnergy( std::sqrt( pp*pp + tarmas*tarmas ) );
00187
00188 theParticleChange.SetEnergyChange( modifiedOriginal.GetKineticEnergy() );
00189 G4DynamicParticle *pd = new G4DynamicParticle;
00190 pd->SetDefinition( targetParticle.GetDefinition() );
00191 pd->SetMomentum( targetParticle.GetMomentum() );
00192 theParticleChange.AddSecondary( pd );
00193 return;
00194 }
00195
00196 G4FastVector<G4ReactionProduct,4> vec;
00197 G4int vecLen = 0;
00198 vec.Initialize( 0 );
00199
00200 G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
00201 G4double massVec[9];
00202 massVec[0] = targetNucleus.AtomicMass( A+1.0, Z );
00203 massVec[1] = theAtomicMass;
00204 massVec[2] = 0.;
00205 if (Z > 1.0) massVec[2] = targetNucleus.AtomicMass(A, Z-1.0);
00206 massVec[3] = 0.;
00207 if (Z > 1.0 && A > 1.0) massVec[3] = targetNucleus.AtomicMass(A-1.0, Z-1.0 );
00208 massVec[4] = 0.;
00209 if (Z > 1.0 && A > 2.0 && A-2.0 > Z-1.0)
00210 massVec[4] = targetNucleus.AtomicMass( A-2.0, Z-1.0 );
00211 massVec[5] = 0.;
00212 if (Z > 2.0 && A > 3.0 && A-3.0 > Z-2.0)
00213 massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-2.0 );
00214 massVec[6] = 0.;
00215 if (A > 1.0 && A-1.0 > Z) massVec[6] = targetNucleus.AtomicMass(A-1.0, Z);
00216 massVec[7] = massVec[3];
00217 massVec[8] = 0.;
00218 if (Z > 2.0 && A > 1.0) massVec[8] = targetNucleus.AtomicMass( A-1.0,Z-2.0 );
00219
00220 twoBody.NuclearReaction(vec, vecLen, originalIncident,
00221 targetNucleus, theAtomicMass, massVec );
00222
00223 theParticleChange.SetStatusChange( stopAndKill );
00224 theParticleChange.SetEnergyChange( 0.0 );
00225
00226 G4DynamicParticle* pd;
00227 for( G4int i=0; i<vecLen; ++i ) {
00228 pd = new G4DynamicParticle();
00229 pd->SetDefinition( vec[i]->GetDefinition() );
00230 pd->SetMomentum( vec[i]->GetMomentum() );
00231 theParticleChange.AddSecondary( pd );
00232 delete vec[i];
00233 }
00234 }
00235
00236
00237
00238
00239
00240
00241
00242
00243 void
00244 G4RPGNeutronInelastic::InitialCollision(G4FastVector<G4ReactionProduct,256>& vec,
00245 G4int& vecLen,
00246 G4ReactionProduct& currentParticle,
00247 G4ReactionProduct& targetParticle,
00248 G4bool& incidentHasChanged,
00249 G4bool& targetHasChanged)
00250 {
00251 G4double KE = currentParticle.GetKineticEnergy()/GeV;
00252
00253 G4int mult;
00254 G4int partType;
00255 std::vector<G4int> fsTypes;
00256 G4int part1;
00257 G4int part2;
00258
00259 G4double testCharge;
00260 G4double testBaryon;
00261 G4double testStrange;
00262
00263
00264
00265 if (targetParticle.GetDefinition() == particleDef[neu]) {
00266 mult = GetMultiplicityT1(KE);
00267 fsTypes = GetFSPartTypesForNN(mult, KE);
00268
00269 part1 = fsTypes[0];
00270 part2 = fsTypes[1];
00271 currentParticle.SetDefinition(particleDef[part1]);
00272 targetParticle.SetDefinition(particleDef[part2]);
00273 if (part1 == pro) {
00274 if (part2 == neu) {
00275 if (G4UniformRand() > 0.5) {
00276 incidentHasChanged = true;
00277 } else {
00278 targetHasChanged = true;
00279 currentParticle.SetDefinition(particleDef[part2]);
00280 targetParticle.SetDefinition(particleDef[part1]);
00281 }
00282 } else {
00283 targetHasChanged = true;
00284 incidentHasChanged = true;
00285 }
00286
00287 } else {
00288 if (part2 > neu && part2 < xi0) targetHasChanged = true;
00289 }
00290
00291 testCharge = 0.0;
00292 testBaryon = 2.0;
00293 testStrange = 0.0;
00294
00295 } else {
00296 mult = GetMultiplicityT0(KE);
00297 fsTypes = GetFSPartTypesForNP(mult, KE);
00298
00299 part1 = fsTypes[0];
00300 part2 = fsTypes[1];
00301 currentParticle.SetDefinition(particleDef[part1]);
00302 targetParticle.SetDefinition(particleDef[part2]);
00303 if (part1 == pro) {
00304 if (part2 == pro) {
00305 incidentHasChanged = true;
00306 } else if (part2 == neu) {
00307 if (G4UniformRand() > 0.5) {
00308 incidentHasChanged = true;
00309 targetHasChanged = true;
00310 } else {
00311 currentParticle.SetDefinition(particleDef[part2]);
00312 targetParticle.SetDefinition(particleDef[part1]);
00313 }
00314
00315 } else if (part2 > neu && part2 < xi0) {
00316 incidentHasChanged = true;
00317 targetHasChanged = true;
00318 }
00319
00320 } else {
00321 targetHasChanged = true;
00322 }
00323
00324 testCharge = 1.0;
00325 testBaryon = 2.0;
00326 testStrange = 0.0;
00327 }
00328
00329
00330
00331
00332
00333
00334 fsTypes.erase(fsTypes.begin());
00335 fsTypes.erase(fsTypes.begin());
00336
00337
00338
00339 G4ReactionProduct* rp(0);
00340 for(G4int i=0; i < mult-2; ++i ) {
00341 partType = fsTypes[i];
00342 rp = new G4ReactionProduct();
00343 rp->SetDefinition(particleDef[partType]);
00344 (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
00345 vec.SetElement(vecLen++, rp);
00346 }
00347
00348
00349
00350 CheckQnums(vec, vecLen, currentParticle, targetParticle,
00351 testCharge, testBaryon, testStrange);
00352
00353 return;
00354 }