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 "G4RPGPiMinusInelastic.hh"
00030 #include "G4SystemOfUnits.hh"
00031 #include "Randomize.hh"
00032
00033 G4HadFinalState*
00034 G4RPGPiMinusInelastic::ApplyYourself(const G4HadProjectile& aTrack,
00035 G4Nucleus& targetNucleus)
00036 {
00037 const G4HadProjectile* originalIncident = &aTrack;
00038
00039 if (originalIncident->GetKineticEnergy()<= 0.1) {
00040 theParticleChange.SetStatusChange(isAlive);
00041 theParticleChange.SetEnergyChange(aTrack.GetKineticEnergy());
00042 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00043 return &theParticleChange;
00044 }
00045
00046
00047
00048 G4DynamicParticle* originalTarget = targetNucleus.ReturnTargetParticle();
00049 G4ReactionProduct targetParticle( originalTarget->GetDefinition() );
00050
00051 G4ReactionProduct currentParticle(
00052 const_cast<G4ParticleDefinition *>(originalIncident->GetDefinition() ) );
00053 currentParticle.SetMomentum( originalIncident->Get4Momentum().vect() );
00054 currentParticle.SetKineticEnergy( originalIncident->GetKineticEnergy() );
00055
00056
00057
00058
00059 G4double ek = originalIncident->GetKineticEnergy();
00060 G4double amas = originalIncident->GetDefinition()->GetPDGMass();
00061
00062 G4double tkin = targetNucleus.Cinema( ek );
00063 ek += tkin;
00064 currentParticle.SetKineticEnergy( ek );
00065 G4double et = ek + amas;
00066 G4double p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00067 G4double pp = currentParticle.GetMomentum().mag();
00068 if( pp > 0.0 ) {
00069 G4ThreeVector momentum = currentParticle.GetMomentum();
00070 currentParticle.SetMomentum( momentum * (p/pp) );
00071 }
00072
00073
00074
00075 tkin = targetNucleus.EvaporationEffects( ek );
00076 ek -= tkin;
00077 currentParticle.SetKineticEnergy( ek );
00078 et = ek + amas;
00079 p = std::sqrt( std::abs((et-amas)*(et+amas)) );
00080 pp = currentParticle.GetMomentum().mag();
00081 if( pp > 0.0 ) {
00082 G4ThreeVector momentum = currentParticle.GetMomentum();
00083 currentParticle.SetMomentum( momentum * (p/pp) );
00084 }
00085
00086 G4ReactionProduct modifiedOriginal = currentParticle;
00087
00088 currentParticle.SetSide( 1 );
00089 targetParticle.SetSide( -1 );
00090 G4bool incidentHasChanged = false;
00091 G4bool targetHasChanged = false;
00092 G4bool quasiElastic = false;
00093 G4FastVector<G4ReactionProduct,256> vec;
00094 G4int vecLen = 0;
00095 vec.Initialize( 0 );
00096
00097 const G4double cutOff = 0.1;
00098 if( currentParticle.GetKineticEnergy() > cutOff )
00099 InitialCollision(vec, vecLen, currentParticle, targetParticle,
00100 incidentHasChanged, targetHasChanged);
00101
00102 CalculateMomenta(vec, vecLen,
00103 originalIncident, originalTarget, modifiedOriginal,
00104 targetNucleus, currentParticle, targetParticle,
00105 incidentHasChanged, targetHasChanged, quasiElastic);
00106
00107 SetUpChange(vec, vecLen,
00108 currentParticle, targetParticle,
00109 incidentHasChanged);
00110
00111 delete originalTarget;
00112 return &theParticleChange;
00113 }
00114
00115
00116
00117
00118
00119
00120
00121
00122 void
00123 G4RPGPiMinusInelastic::InitialCollision(G4FastVector<G4ReactionProduct,256>& vec,
00124 G4int& vecLen,
00125 G4ReactionProduct& currentParticle,
00126 G4ReactionProduct& targetParticle,
00127 G4bool& incidentHasChanged,
00128 G4bool& targetHasChanged)
00129 {
00130 G4double KE = currentParticle.GetKineticEnergy()/GeV;
00131
00132 G4int mult;
00133 G4int partType;
00134 std::vector<G4int> fsTypes;
00135
00136 G4double testCharge;
00137 G4double testBaryon;
00138 G4double testStrange;
00139
00140
00141
00142 if (targetParticle.GetDefinition() == particleDef[pro]) {
00143 mult = GetMultiplicityT12(KE);
00144 fsTypes = GetFSPartTypesForPimP(mult, KE);
00145 partType = fsTypes[0];
00146 if (partType != pro) {
00147 targetHasChanged = true;
00148 targetParticle.SetDefinition(particleDef[partType]);
00149 }
00150
00151 testCharge = 0.0;
00152 testBaryon = 1.0;
00153 testStrange = 0.0;
00154
00155 } else {
00156 mult = GetMultiplicityT32(KE);
00157 fsTypes = GetFSPartTypesForPimN(mult, KE);
00158 partType = fsTypes[0];
00159 if (partType != neu) {
00160 targetHasChanged = true;
00161 targetParticle.SetDefinition(particleDef[partType]);
00162 }
00163
00164 testCharge = -1.0;
00165 testBaryon = 1.0;
00166 testStrange = 0.0;
00167 }
00168
00169
00170
00171 fsTypes.erase(fsTypes.begin());
00172
00173
00174
00175 G4int choose = -1;
00176 for(G4int i=0; i < mult-1; ++i ) {
00177 partType = fsTypes[i];
00178 if (partType == pim) {
00179 choose = i;
00180 break;
00181 }
00182 }
00183 if (choose == -1) {
00184 incidentHasChanged = true;
00185 choose = G4int(G4UniformRand()*(mult-1) );
00186 partType = fsTypes[choose];
00187 currentParticle.SetDefinition(particleDef[partType]);
00188 }
00189
00190 fsTypes.erase(fsTypes.begin()+choose);
00191
00192
00193
00194 G4ReactionProduct* rp(0);
00195 for(G4int i=0; i < mult-2; ++i ) {
00196 partType = fsTypes[i];
00197 rp = new G4ReactionProduct();
00198 rp->SetDefinition(particleDef[partType]);
00199 (G4UniformRand() < 0.5) ? rp->SetSide(-1) : rp->SetSide(1);
00200 if (partType > pim && partType < pro) rp->SetMayBeKilled(false);
00201 vec.SetElement(vecLen++, rp);
00202 }
00203
00204
00205
00206
00207
00208
00209 CheckQnums(vec, vecLen, currentParticle, targetParticle,
00210 testCharge, testBaryon, testStrange);
00211
00212 return;
00213 }