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 "G4RPGProtonInelastic.hh"
00030 #include "G4SystemOfUnits.hh"
00031 #include "Randomize.hh"
00032
00033 G4HadFinalState*
00034 G4RPGProtonInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00035 G4Nucleus& targetNucleus )
00036 {
00037 theParticleChange.Clear();
00038 const G4HadProjectile *originalIncident = &aTrack;
00039 if (originalIncident->GetKineticEnergy()<= 0.1)
00040 {
00041 theParticleChange.SetStatusChange(isAlive);
00042 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00043 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00044 return &theParticleChange;
00045 }
00046
00047
00048
00049
00050 G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
00051
00052 if (originalIncident->GetKineticEnergy()/GeV < 0.01+2.*G4UniformRand()/9. )
00053 {
00054 SlowProton( originalIncident, targetNucleus );
00055 delete originalTarget;
00056 return &theParticleChange;
00057 }
00058
00059
00060
00061
00062 G4double ek = originalIncident->GetKineticEnergy();
00063 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00064 G4ReactionProduct modifiedOriginal;
00065 modifiedOriginal = *originalIncident;
00066
00067 G4double tkin = targetNucleus.Cinema( ek );
00068 ek += tkin;
00069 modifiedOriginal.SetKineticEnergy(ek);
00070 G4double et = ek + amas;
00071 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00072 G4double pp = modifiedOriginal.GetMomentum().mag();
00073 if (pp > 0.0) {
00074 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00075 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00076 }
00077
00078
00079
00080 tkin = targetNucleus.EvaporationEffects(ek);
00081 ek -= tkin;
00082 modifiedOriginal.SetKineticEnergy(ek);
00083 et = ek + amas;
00084 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00085 pp = modifiedOriginal.GetMomentum().mag();
00086 if (pp > 0.0) {
00087 G4ThreeVector momentum = modifiedOriginal.GetMomentum();
00088 modifiedOriginal.SetMomentum( momentum * (p/pp) );
00089 }
00090 const G4double cutOff = 0.1;
00091 if (modifiedOriginal.GetKineticEnergy() < cutOff) {
00092 SlowProton( originalIncident, targetNucleus );
00093 delete originalTarget;
00094 return &theParticleChange;
00095 }
00096
00097 G4ReactionProduct currentParticle = modifiedOriginal;
00098 G4ReactionProduct targetParticle;
00099 targetParticle = *originalTarget;
00100 currentParticle.SetSide( 1 );
00101 targetParticle.SetSide( -1 );
00102 G4bool incidentHasChanged = false;
00103 G4bool targetHasChanged = false;
00104 G4bool quasiElastic = false;
00105 G4FastVector<G4ReactionProduct,256> vec;
00106 G4int vecLen = 0;
00107 vec.Initialize( 0 );
00108
00109 InitialCollision(vec, vecLen, currentParticle, targetParticle,
00110 incidentHasChanged, targetHasChanged);
00111
00112 CalculateMomenta(vec, vecLen,
00113 originalIncident, originalTarget, modifiedOriginal,
00114 targetNucleus, currentParticle, targetParticle,
00115 incidentHasChanged, targetHasChanged, quasiElastic);
00116
00117 SetUpChange( vec, vecLen,
00118 currentParticle, targetParticle,
00119 incidentHasChanged );
00120
00121 delete originalTarget;
00122 return &theParticleChange;
00123 }
00124
00125
00126 void
00127 G4RPGProtonInelastic::SlowProton(const G4HadProjectile *originalIncident,
00128 G4Nucleus &targetNucleus )
00129 {
00130 const G4double A = targetNucleus.GetA_asInt();
00131 const G4double Z = targetNucleus.GetZ_asInt();
00132
00133
00134
00135 G4double theAtomicMass = targetNucleus.AtomicMass( A, Z );
00136 G4double massVec[9];
00137 massVec[0] = targetNucleus.AtomicMass( A+1.0, Z+1.0 );
00138 massVec[1] = 0.;
00139 if (A > Z+1.0)
00140 massVec[1] = targetNucleus.AtomicMass( A , Z+1.0 );
00141 massVec[2] = theAtomicMass;
00142 massVec[3] = 0.;
00143 if (A > 1.0 && A-1.0 > Z)
00144 massVec[3] = targetNucleus.AtomicMass( A-1.0, Z );
00145 massVec[4] = 0.;
00146 if (A > 2.0 && A-2.0 > Z)
00147 massVec[4] = targetNucleus.AtomicMass( A-2.0, Z );
00148 massVec[5] = 0.;
00149 if (A > 3.0 && Z > 1.0 && A-3.0 > Z-1.0)
00150 massVec[5] = targetNucleus.AtomicMass( A-3.0, Z-1.0 );
00151 massVec[6] = 0.;
00152 if (A > 1.0 && A-1.0 > Z+1.0)
00153 massVec[6] = targetNucleus.AtomicMass( A-1.0, Z+1.0 );
00154 massVec[7] = massVec[3];
00155 massVec[8] = 0.;
00156 if (A > 1.0 && Z > 1.0)
00157 massVec[8] = targetNucleus.AtomicMass( A-1.0, Z-1.0 );
00158
00159 G4FastVector<G4ReactionProduct,4> vec;
00160 G4int vecLen = 0;
00161 vec.Initialize( 0 );
00162
00163 twoBody.NuclearReaction( vec, vecLen, originalIncident,
00164 targetNucleus, theAtomicMass, massVec );
00165
00166 theParticleChange.SetStatusChange( stopAndKill );
00167 theParticleChange.SetEnergyChange( 0.0 );
00168
00169 G4DynamicParticle *pd;
00170 for( G4int i=0; i<vecLen; ++i )
00171 {
00172 pd = new G4DynamicParticle();
00173 pd->SetDefinition( vec[i]->GetDefinition() );
00174 pd->SetMomentum( vec[i]->GetMomentum() );
00175 theParticleChange.AddSecondary( pd );
00176 delete vec[i];
00177 }
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187 void
00188 G4RPGProtonInelastic::InitialCollision(G4FastVector<G4ReactionProduct,256>& vec,
00189 G4int& vecLen,
00190 G4ReactionProduct& currentParticle,
00191 G4ReactionProduct& targetParticle,
00192 G4bool& incidentHasChanged,
00193 G4bool& targetHasChanged)
00194 {
00195 G4double KE = currentParticle.GetKineticEnergy()/GeV;
00196
00197 G4int mult;
00198 G4int partType;
00199 std::vector<G4int> fsTypes;
00200 G4int part1;
00201 G4int part2;
00202
00203 G4double testCharge;
00204 G4double testBaryon;
00205 G4double testStrange;
00206
00207
00208
00209 if (targetParticle.GetDefinition() == particleDef[pro]) {
00210 mult = GetMultiplicityT1(KE);
00211 fsTypes = GetFSPartTypesForPP(mult, KE);
00212
00213 part1 = fsTypes[0];
00214 part2 = fsTypes[1];
00215 currentParticle.SetDefinition(particleDef[part1]);
00216 targetParticle.SetDefinition(particleDef[part2]);
00217 if (part1 == pro) {
00218 if (part2 == neu) {
00219 if (G4UniformRand() > 0.5) {
00220 incidentHasChanged = true;
00221 targetParticle.SetDefinition(particleDef[part1]);
00222 currentParticle.SetDefinition(particleDef[part2]);
00223 } else {
00224 targetHasChanged = true;
00225 }
00226 } else if (part2 > neu && part2 < xi0) {
00227 targetHasChanged = true;
00228 }
00229
00230 } else {
00231 targetHasChanged = true;
00232 incidentHasChanged = true;
00233 }
00234
00235 testCharge = 2.0;
00236 testBaryon = 2.0;
00237 testStrange = 0.0;
00238
00239 } else {
00240 mult = GetMultiplicityT0(KE);
00241 fsTypes = GetFSPartTypesForPN(mult, KE);
00242
00243 part1 = fsTypes[0];
00244 part2 = fsTypes[1];
00245 currentParticle.SetDefinition(particleDef[part1]);
00246 targetParticle.SetDefinition(particleDef[part2]);
00247 if (part1 == pro) {
00248 if (part2 == pro) {
00249 targetHasChanged = true;
00250 } else if (part2 == neu) {
00251 if (G4UniformRand() > 0.5) {
00252 incidentHasChanged = true;
00253 targetHasChanged = true;
00254 targetParticle.SetDefinition(particleDef[part1]);
00255 currentParticle.SetDefinition(particleDef[part2]);
00256 }
00257 } else {
00258 targetHasChanged = true;
00259 }
00260
00261 } else {
00262 incidentHasChanged = true;
00263 if (part2 > neu && part2 < xi0) targetHasChanged = true;
00264 }
00265
00266 testCharge = 1.0;
00267 testBaryon = 2.0;
00268 testStrange = 0.0;
00269 }
00270
00271
00272
00273 fsTypes.erase(fsTypes.begin());
00274 fsTypes.erase(fsTypes.begin());
00275
00276
00277
00278 G4ReactionProduct* rp(0);
00279 for(G4int i=0; i < mult-2; ++i ) {
00280 partType = fsTypes[i];
00281 rp = new G4ReactionProduct();
00282 rp->SetDefinition(particleDef[partType]);
00283 (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
00284 vec.SetElement(vecLen++, rp);
00285 }
00286
00287
00288
00289 CheckQnums(vec, vecLen, currentParticle, targetParticle,
00290 testCharge, testBaryon, testStrange);
00291
00292 return;
00293 }