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
00040
00041
00042 #include "G4ChargeExchangeProcess.hh"
00043 #include "globals.hh"
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4CrossSectionDataStore.hh"
00046 #include "G4HadronElasticDataSet.hh"
00047 #include "G4Element.hh"
00048 #include "G4ElementVector.hh"
00049 #include "G4IsotopeVector.hh"
00050 #include "G4Neutron.hh"
00051 #include "G4Proton.hh"
00052 #include "G4PhysicsLinearVector.hh"
00053
00054 G4ChargeExchangeProcess::G4ChargeExchangeProcess(const G4String& procName)
00055 : G4HadronicProcess(procName,fChargeExchange), first(true)
00056 {
00057 thEnergy = 20.*MeV;
00058 pPDG = 0;
00059 verboseLevel= 1;
00060 AddDataSet(new G4HadronElasticDataSet);
00061 theProton = G4Proton::Proton();
00062 theNeutron = G4Neutron::Neutron();
00063 theAProton = G4AntiProton::AntiProton();
00064 theANeutron = G4AntiNeutron::AntiNeutron();
00065 thePiPlus = G4PionPlus::PionPlus();
00066 thePiMinus = G4PionMinus::PionMinus();
00067 thePiZero = G4PionZero::PionZero();
00068 theKPlus = G4KaonPlus::KaonPlus();
00069 theKMinus = G4KaonMinus::KaonMinus();
00070 theK0S = G4KaonZeroShort::KaonZeroShort();
00071 theK0L = G4KaonZeroLong::KaonZeroLong();
00072 theL = G4Lambda::Lambda();
00073 theAntiL = G4AntiLambda::AntiLambda();
00074 theSPlus = G4SigmaPlus::SigmaPlus();
00075 theASPlus = G4AntiSigmaPlus::AntiSigmaPlus();
00076 theSMinus = G4SigmaMinus::SigmaMinus();
00077 theASMinus = G4AntiSigmaMinus::AntiSigmaMinus();
00078 theS0 = G4SigmaZero::SigmaZero();
00079 theAS0 = G4AntiSigmaZero::AntiSigmaZero();
00080 theXiMinus = G4XiMinus::XiMinus();
00081 theXi0 = G4XiZero::XiZero();
00082 theAXiMinus = G4AntiXiMinus::AntiXiMinus();
00083 theAXi0 = G4AntiXiZero::AntiXiZero();
00084 theOmega = G4OmegaMinus::OmegaMinus();
00085 theAOmega = G4AntiOmegaMinus::AntiOmegaMinus();
00086 theD = G4Deuteron::Deuteron();
00087 theT = G4Triton::Triton();
00088 theA = G4Alpha::Alpha();
00089 theHe3 = G4He3::He3();
00090 }
00091
00092 G4ChargeExchangeProcess::~G4ChargeExchangeProcess()
00093 {
00094 delete factors;
00095 }
00096
00097 void G4ChargeExchangeProcess::
00098 BuildPhysicsTable(const G4ParticleDefinition& aParticleType)
00099 {
00100 if(first) {
00101 first = false;
00102 theParticle = &aParticleType;
00103 pPDG = theParticle->GetPDGEncoding();
00104
00105 store = G4HadronicProcess::GetCrossSectionDataStore();
00106
00107 const size_t n = 10;
00108 if(theParticle == thePiPlus || theParticle == thePiMinus ||
00109 theParticle == theKPlus || theParticle == theKMinus ||
00110 theParticle == theK0S || theParticle == theK0L) {
00111
00112 G4double F[n] = {0.33,0.27,0.29,0.31,0.27,0.18,0.13,0.1,0.09,0.07};
00113 factors = new G4PhysicsLinearVector(0.0,2.0*GeV,n);
00114 for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
00115
00116 } else {
00117
00118 G4double F[n] = {0.50,0.45,0.40,0.35,0.30,0.25,0.06,0.04,0.005,0.0};
00119 factors = new G4PhysicsLinearVector(0.0,4.0*GeV,n);
00120 for(size_t i=0; i<n; i++) {factors->PutValue(i,F[i]);}
00121 }
00122
00123
00124 if(verboseLevel>1)
00125 G4cout << "G4ChargeExchangeProcess for "
00126 << theParticle->GetParticleName()
00127 << G4endl;
00128 }
00129 G4HadronicProcess::BuildPhysicsTable(aParticleType);
00130 }
00131
00132 G4double G4ChargeExchangeProcess::GetElementCrossSection(
00133 const G4DynamicParticle* dp,
00134 const G4Element* elm,
00135 const G4Material* mat)
00136 {
00137
00138 G4double Z = elm->GetZ();
00139 G4int iz = G4int(Z);
00140 G4double x = 0.0;
00141
00142
00143 if(iz == 1 || dp->GetKineticEnergy() < thEnergy) return x;
00144
00145 if(verboseLevel>1)
00146 G4cout << "G4ChargeExchangeProcess compute GHAD CS for element "
00147 << elm->GetName()
00148 << G4endl;
00149 x = store->GetCrossSection(dp, elm, mat);
00150
00151 if(verboseLevel>1)
00152 G4cout << "G4ChargeExchangeProcess cross(mb)= " << x/millibarn
00153 << " E(MeV)= " << dp->GetKineticEnergy()
00154 << " " << theParticle->GetParticleName()
00155 << " in Z= " << iz
00156 << G4endl;
00157 G4bool b;
00158 G4double A = elm->GetN();
00159 G4double ptot = dp->GetTotalMomentum();
00160 x *= factors->GetValue(ptot, b)/std::pow(A, 0.42);
00161 if(theParticle == thePiPlus || theParticle == theProton ||
00162 theParticle == theKPlus || theParticle == theANeutron)
00163 { x *= (1.0 - Z/A); }
00164
00165 else if(theParticle == thePiMinus || theParticle == theNeutron ||
00166 theParticle == theKMinus || theParticle == theAProton)
00167 { x *= Z/A; }
00168
00169 if(theParticle->GetPDGMass() < GeV) {
00170 if(ptot > 2.*GeV) x *= 4.0*GeV*GeV/(ptot*ptot);
00171 }
00172
00173 if(verboseLevel>1)
00174 G4cout << "Corrected cross(mb)= " << x/millibarn << G4endl;
00175
00176 return x;
00177 }
00178
00179 G4bool G4ChargeExchangeProcess::
00180 IsApplicable(const G4ParticleDefinition& aParticleType)
00181 {
00182 const G4ParticleDefinition* p = &aParticleType;
00183 return (p == thePiPlus || p == thePiMinus ||
00184 p == theProton || p == theNeutron ||
00185 p == theAProton|| p == theANeutron||
00186 p == theKPlus || p == theKMinus ||
00187 p == theK0S || p == theK0L ||
00188 p == theL);
00189 }
00190
00191 void G4ChargeExchangeProcess::
00192 DumpPhysicsTable(const G4ParticleDefinition& aParticleType)
00193 {
00194 store->DumpPhysicsTable(aParticleType);
00195 }