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
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #include "G4ChargeExchange.hh"
00040 #include "G4PhysicalConstants.hh"
00041 #include "G4SystemOfUnits.hh"
00042 #include "G4ParticleTable.hh"
00043 #include "G4ParticleDefinition.hh"
00044 #include "G4IonTable.hh"
00045 #include "Randomize.hh"
00046 #include "G4NucleiProperties.hh"
00047
00048 G4ChargeExchange::G4ChargeExchange() : G4HadronicInteraction("Charge Exchange")
00049 {
00050 SetMinEnergy( 0.0*GeV );
00051 SetMaxEnergy( 100.*TeV );
00052
00053 lowEnergyRecoilLimit = 100.*keV;
00054 lowestEnergyLimit = 1.*MeV;
00055
00056 theProton = G4Proton::Proton();
00057 theNeutron = G4Neutron::Neutron();
00058 theAProton = G4AntiProton::AntiProton();
00059 theANeutron = G4AntiNeutron::AntiNeutron();
00060 thePiPlus = G4PionPlus::PionPlus();
00061 thePiMinus = G4PionMinus::PionMinus();
00062 thePiZero = G4PionZero::PionZero();
00063 theKPlus = G4KaonPlus::KaonPlus();
00064 theKMinus = G4KaonMinus::KaonMinus();
00065 theK0S = G4KaonZeroShort::KaonZeroShort();
00066 theK0L = G4KaonZeroLong::KaonZeroLong();
00067 theL = G4Lambda::Lambda();
00068 theAntiL = G4AntiLambda::AntiLambda();
00069 theSPlus = G4SigmaPlus::SigmaPlus();
00070 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
00071 theSMinus = G4SigmaMinus::SigmaMinus();
00072 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
00073 theS0 = G4SigmaZero::SigmaZero();
00074 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
00075 theXiMinus = G4XiMinus::XiMinus();
00076 theXi0 = G4XiZero::XiZero();
00077 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00078 theAXi0 = G4AntiXiZero::AntiXiZero();
00079 theOmega = G4OmegaMinus::OmegaMinus();
00080 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
00081 theD = G4Deuteron::Deuteron();
00082 theT = G4Triton::Triton();
00083 theA = G4Alpha::Alpha();
00084 theHe3 = G4He3::He3();
00085 }
00086
00087 G4ChargeExchange::~G4ChargeExchange()
00088 {}
00089
00090 G4HadFinalState* G4ChargeExchange::ApplyYourself(
00091 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
00092 {
00093 theParticleChange.Clear();
00094 const G4HadProjectile* aParticle = &aTrack;
00095 G4double ekin = aParticle->GetKineticEnergy();
00096
00097 G4int A = targetNucleus.GetA_asInt();
00098 G4int Z = targetNucleus.GetZ_asInt();
00099
00100 if(ekin <= lowestEnergyLimit || A < 3) {
00101 theParticleChange.SetEnergyChange(ekin);
00102 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00103 return &theParticleChange;
00104 }
00105
00106 G4double plab = aParticle->GetTotalMomentum();
00107
00108 if (verboseLevel > 1)
00109 G4cout << "G4ChargeExchange::DoIt: Incident particle plab="
00110 << plab/GeV << " GeV/c "
00111 << " ekin(MeV) = " << ekin/MeV << " "
00112 << aParticle->GetDefinition()->GetParticleName() << G4endl;
00113
00114
00115 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
00116
00117 G4int N = A - Z;
00118 G4int projPDG = theParticle->GetPDGEncoding();
00119 if (verboseLevel > 1)
00120 G4cout << "G4ChargeExchange for " << theParticle->GetParticleName()
00121 << " PDGcode= " << projPDG << " on nucleus Z= " << Z
00122 << " A= " << A << " N= " << N
00123 << G4endl;
00124
00125 G4ParticleDefinition * theDef = 0;
00126
00127 G4double mass2 = G4NucleiProperties::GetNuclearMass((G4double)A, (G4double)Z);
00128 G4LorentzVector lv1 = aParticle->Get4Momentum();
00129 G4LorentzVector lv0(0.0,0.0,0.0,mass2);
00130
00131 G4LorentzVector lv = lv0 + lv1;
00132 G4ThreeVector bst = lv.boostVector();
00133 lv1.boost(-bst);
00134 lv0.boost(-bst);
00135
00136
00137 G4bool theHyperon = false;
00138 G4ParticleDefinition* theRecoil = 0;
00139 G4ParticleDefinition* theSecondary = 0;
00140
00141 if(theParticle == theProton) {
00142 theSecondary = theNeutron;
00143 Z++;
00144 } else if(theParticle == theNeutron) {
00145 theSecondary = theProton;
00146 Z--;
00147 } else if(theParticle == thePiPlus) {
00148 theSecondary = thePiZero;
00149 Z++;
00150 } else if(theParticle == thePiMinus) {
00151 theSecondary = thePiZero;
00152 Z--;
00153 } else if(theParticle == theKPlus) {
00154 if(G4UniformRand()<0.5) theSecondary = theK0S;
00155 else theSecondary = theK0L;
00156 Z++;
00157 } else if(theParticle == theKMinus) {
00158 if(G4UniformRand()<0.5) theSecondary = theK0S;
00159 else theSecondary = theK0L;
00160 Z--;
00161 } else if(theParticle == theK0S || theParticle == theK0L) {
00162 if(G4UniformRand()*A < G4double(Z)) {
00163 theSecondary = theKPlus;
00164 Z--;
00165 } else {
00166 theSecondary = theKMinus;
00167 Z++;
00168 }
00169 } else if(theParticle == theANeutron) {
00170 theSecondary = theAProton;
00171 Z++;
00172 } else if(theParticle == theAProton) {
00173 theSecondary = theANeutron;
00174 Z--;
00175 } else if(theParticle == theL) {
00176 G4double x = G4UniformRand();
00177 if(G4UniformRand()*A < G4double(Z)) {
00178 if(x < 0.2) {
00179 theSecondary = theS0;
00180 } else if (x < 0.4) {
00181 theSecondary = theSPlus;
00182 Z--;
00183 } else if (x < 0.6) {
00184 theSecondary = theProton;
00185 theRecoil = theL;
00186 theHyperon = true;
00187 A--;
00188 } else if (x < 0.8) {
00189 theSecondary = theProton;
00190 theRecoil = theS0;
00191 theHyperon = true;
00192 A--;
00193 } else {
00194 theSecondary = theNeutron;
00195 theRecoil = theSPlus;
00196 theHyperon = true;
00197 A--;
00198 }
00199 } else {
00200 if(x < 0.2) {
00201 theSecondary = theS0;
00202 } else if (x < 0.4) {
00203 theSecondary = theSMinus;
00204 Z++;
00205 } else if (x < 0.6) {
00206 theSecondary = theNeutron;
00207 theRecoil = theL;
00208 A--;
00209 theHyperon = true;
00210 } else if (x < 0.8) {
00211 theSecondary = theNeutron;
00212 theRecoil = theS0;
00213 theHyperon = true;
00214 A--;
00215 } else {
00216 theSecondary = theProton;
00217 theRecoil = theSMinus;
00218 theHyperon = true;
00219 A--;
00220 }
00221 }
00222 }
00223
00224 if (Z == 1 && A == 2) theDef = theD;
00225 else if (Z == 1 && A == 3) theDef = theT;
00226 else if (Z == 2 && A == 3) theDef = theHe3;
00227 else if (Z == 2 && A == 4) theDef = theA;
00228 else {
00229 theDef =
00230 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(Z,A,0.0);
00231 }
00232 if(!theSecondary) { return &theParticleChange; }
00233
00234 G4double m11 = theSecondary->GetPDGMass();
00235 G4double m21 = theDef->GetPDGMass();
00236 if(theRecoil) { m21 += theRecoil->GetPDGMass(); }
00237 else { theRecoil = theDef; }
00238
00239 G4double etot = lv0.e() + lv1.e();
00240
00241
00242 if(etot < m11 + m21) {
00243 theParticleChange.SetEnergyChange(ekin);
00244 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
00245 return &theParticleChange;
00246 }
00247
00248 G4ThreeVector p1 = lv1.vect();
00249 G4double e1 = 0.5*etot*(1.0 - (m21*m21 - m11*m11)/(etot*etot));
00250
00251 G4double ptot = std::sqrt(e1*e1 - m11*m11);
00252
00253 G4double tmax = 4.0*ptot*ptot;
00254 G4double g2 = GeV*GeV;
00255
00256 G4double t = g2*SampleT(tmax/g2, A);
00257
00258 if(verboseLevel>1)
00259 G4cout <<"## G4ChargeExchange t= " << t << " tmax= " << tmax
00260 << " ptot= " << ptot << G4endl;
00261
00262
00263 G4double phi = G4UniformRand()*twopi;
00264 G4double cost = 1. - 2.0*t/tmax;
00265 if(std::abs(cost) > 1.0) cost = 1.0;
00266 G4double sint = std::sqrt((1.0-cost)*(1.0+cost));
00267
00268
00269
00270
00271 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
00272 v1 *= ptot;
00273 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),e1);
00274 G4LorentzVector nlv0 = lv0 + lv1 - nlv1;
00275
00276 nlv0.boost(bst);
00277 nlv1.boost(bst);
00278
00279 theParticleChange.SetStatusChange(stopAndKill);
00280 theParticleChange.SetEnergyChange(0.0);
00281 G4DynamicParticle * aSec = new G4DynamicParticle(theSecondary, nlv1);
00282 theParticleChange.AddSecondary(aSec);
00283
00284 G4double erec = nlv0.e() - m21;
00285
00286
00287
00288 if(theHyperon) {
00289 theParticleChange.SetLocalEnergyDeposit(erec);
00290 aSec = new G4DynamicParticle();
00291 aSec->SetDefinition(theRecoil);
00292 aSec->SetKineticEnergy(0.0);
00293 } else if(erec > lowEnergyRecoilLimit) {
00294 aSec = new G4DynamicParticle(theRecoil, nlv0);
00295 theParticleChange.AddSecondary(aSec);
00296 } else {
00297 if(erec < 0.0) erec = 0.0;
00298 theParticleChange.SetLocalEnergyDeposit(erec);
00299 }
00300 return &theParticleChange;
00301 }
00302
00303 G4double G4ChargeExchange::SampleT(G4double tmax, G4double A)
00304 {
00305 G4double aa, bb, cc, dd;
00306 if (A <= 62.) {
00307 aa = std::pow(A, 1.63);
00308 bb = 14.5*std::pow(A, 0.66);
00309 cc = 1.4*std::pow(A, 0.33);
00310 dd = 10.;
00311 } else {
00312 aa = std::pow(A, 1.33);
00313 bb = 60.*std::pow(A, 0.33);
00314 cc = 0.4*std::pow(A, 0.40);
00315 dd = 10.;
00316 }
00317 G4double x1 = (1.0 - std::exp(-tmax*bb))*aa/bb;
00318 G4double x2 = (1.0 - std::exp(-tmax*dd))*cc/dd;
00319
00320 G4double t;
00321 G4double y = bb;
00322 if(G4UniformRand()*(x1 + x2) < x2) y = dd;
00323
00324 do {t = -std::log(G4UniformRand())/y;} while (t > tmax);
00325
00326 return t;
00327 }
00328