Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // ------------------------------------------- 27 // 28 // GEANT4 Class file 29 // 30 // 31 // File name: G4ChargeExchangeXS 32 // 33 34 #include "G4ChargeExchangeXS.hh" 35 #include "G4DynamicParticle.hh" 36 #include "G4ElementTable.hh" 37 #include "G4Material.hh" 38 #include "G4Element.hh" 39 #include "G4IsotopeList.hh" 40 #include "G4HadronicParameters.hh" 41 #include "Randomize.hh" 42 #include "G4SystemOfUnits.hh" 43 #include "G4NucleiProperties.hh" 44 #include "G4Pow.hh" 45 46 #include "G4PionZero.hh" 47 #include "G4Eta.hh" 48 #include "G4KaonZeroLong.hh" 49 #include "G4KaonZeroShort.hh" 50 #include "G4KaonPlus.hh" 51 #include "G4KaonMinus.hh" 52 #include "G4ParticleTable.hh" 53 54 namespace { 55 // V. Lyubovitsky parameterisation 56 const G4double piA[5] = {430., 36., 1.37, 57 const G4double pAP[5] = {1.04, 1.26, 1.35, 58 const G4double pC0[5] = {12.7, 6.0, 6.84, 59 const G4double pC1[5] = {1.57, 1.6, 1.7, 60 const G4double pG0[5] = {2.55, 4.6, 3.7, 61 const G4double pG1[5] = {-0.23, -0.5, 0., 62 63 const G4double beta_prime_pi = 0.0410; 64 } 65 66 G4ChargeExchangeXS::G4ChargeExchangeXS() 67 { 68 if (verboseLevel > 1) { 69 G4cout << "G4ChargeExchangeXS::G4ChargeEx 70 } 71 g4calc = G4Pow::GetInstance(); 72 auto table = G4ParticleTable::GetParticleTab 73 const G4String nam[5] = {"pi0", "eta", "eta_ 74 for (G4int i=0; i<5; ++i) { 75 fPionSecPD[i] = table->FindParticle(nam[i] 76 if (nullptr == fPionSecPD[i]) { 77 G4ExceptionDescription ed; 78 ed << "### meson " << nam[i] << " is not 79 G4Exception("G4ChargeExchangeXS::G4Charg 80 FatalException, ed,""); 81 } 82 } 83 } 84 85 // Print the information of this .cc file 86 void G4ChargeExchangeXS::CrossSectionDescripti 87 { 88 outFile << "G4ChargeExchangeXS calculates ch 89 << "pi+, pi-, K+, K-, KL\n"; 90 } 91 92 G4bool G4ChargeExchangeXS::IsElementApplicable 93 94 { 95 return true; 96 } 97 98 G4double 99 G4ChargeExchangeXS::GetElementCrossSection(con 100 G4int ZZ, const G4Material* 101 { 102 G4double result = 0.0; 103 const G4double pE = aParticle->GetTotalEnerg 104 if (pE <= fEnergyLimit) { return result; } 105 106 auto part = aParticle->GetDefinition(); 107 G4int pdg = part->GetPDGEncoding(); 108 109 // Get or calculate the proton mass, particl 110 G4double tM = CLHEP::proton_mass_c2; 111 G4double pM = part->GetPDGMass(); 112 G4double lorentz_s = tM*tM + 2*tM*pE + pM*p 113 if (lorentz_s <= (tM + pM)*(tM + pM)) { retu 114 115 const G4int Z = std::min(ZZ, ZMAXNUCLEARDATA 116 const G4int A = G4lrint(aeff[Z]); 117 118 if (verboseLevel > 1) { 119 G4cout << "### G4ChargeExchangeXS: " << pa 120 << " Z=" << Z << " A=" << A << " Etot(GeV 121 << " s(GeV^2)=" << lorentz_s/(CLHEP::GeV* 122 } 123 124 // For unit conversion 125 const G4double inv1e7 = 0.1/(CLHEP::GeV*CLHE 126 const G4double fact = 1e-30*CLHEP::cm2; 127 const G4double pfact = 0.1/CLHEP::GeV; 128 const G4double kfact = 56.3*fact; 129 const G4double csmax = 1e-16; 130 131 // The approximation of Glauber-Gribov formu 132 // proton to nuclei Z^(2/3). The factor g4ca 133 // takes into account absorption of pi0 and 134 135 // pi- + p -> n + meson (0- pi0, 1- eta, 2- 136 if (pdg == -211) { 137 G4double z23 = g4calc->Z23(Z); 138 G4double x = lorentz_s*inv1e7; 139 G4double logX = G4Log(x); 140 G4double logA = g4calc->logZ(A); 141 G4double xf = g4calc->powZ(A, -beta_prime_ 142 G4double sum = 0.0; 143 for (G4int i=0; i<5; ++i) { 144 G4double xg = std::max(1.0 + pG0[i] + pG 145 G4double xc = std::max(pC0[i] + pC1[i]*l 146 G4double xs = z23*piA[i]*g4calc->powA(x, 147 sum += xs; 148 fXSecPion[i] = sum; 149 } 150 result = sum*fact; 151 } 152 153 // pi+ + n -> p + meson (0- pi0, 1- eta, 2- 154 else if (pdg == 211) { 155 G4double n23 = g4calc->Z23(A - Z); 156 G4double x = lorentz_s*inv1e7; 157 G4double logX = G4Log(x); 158 G4double logA = g4calc->logZ(A); 159 G4double xf = g4calc->powZ(A, -beta_prime_ 160 161 // hydrogen target case Z = A = 1 162 // the cross section is defined by fractio 163 if (1 == Z) { n23 = ComputeDeuteronFractio 164 G4double sum = 0.0; 165 for (G4int i=0; i<5; ++i) { 166 G4double xg = std::max(1.0 + pG0[i] + pG 167 G4double xc = std::max(pC0[i] + pC1[i]*l 168 G4double xs = n23*piA[i]*g4calc->powA(x, 169 sum += xs; 170 fXSecPion[i] = sum; 171 } 172 result = sum*fact; 173 } 174 175 // Kaon x-sections depend on the primary par 176 // K- + p -> Kbar + n 177 else if (pdg == -321) { 178 G4double p_momentum = std::sqrt(pE*pE - pM 179 result = g4calc->Z23(Z)*g4calc->powA(p_mom 180 } 181 182 // K+ + n -> Kbar + p 183 else if (pdg == 321) { 184 G4double p_momentum = std::sqrt(pE*pE - pM 185 G4double n23 = g4calc->Z23(A-Z); 186 // hydrogen target case Z = A = 1 187 // the cross section is defined by fractio 188 if (1 == Z) { n23 = ComputeDeuteronFractio 189 result = n23*g4calc->powA(p_momentum, -1.6 190 } 191 192 // KL 193 else if (pdg == 130) { 194 // Cross section of KL = 0.5*(Cross sectio 195 const G4double p_momentum = std::sqrt(pE*p 196 result = 0.5*(g4calc->Z23(Z) + g4calc->Z23 197 g4calc->powA(p_momentum, -1.60)*kfact; 198 } 199 result *= fFactor; 200 if (verboseLevel > 1) { 201 G4cout << " Done for " << part->GetPart 202 << " res(mb)=" << result/CLHEP::millibar 203 } 204 return result; 205 } 206 207 const G4ParticleDefinition* 208 G4ChargeExchangeXS::SampleSecondaryType(const 209 const 210 { 211 const G4ParticleDefinition* pd = nullptr; 212 G4int pdg = std::abs(part->GetPDGEncoding()) 213 214 // pi- + p / pi+ + n 215 if (pdg == 211) { 216 pd = fPionSecPD[0]; 217 G4double x = fXSecPion[4]*G4UniformRand(); 218 for (G4int i=0; i<5; ++i) { 219 if (x <= fXSecPion[i]) { 220 pd = fPionSecPD[i]; 221 break; 222 } 223 } 224 } 225 226 // K- + p / K+ + n 227 // Equal opportunity of producing k-short an 228 else if (pdg == 321) { 229 if (G4UniformRand() >= 0.5) { 230 pd = G4KaonZeroLong::KaonZeroLong(); 231 } 232 else { 233 pd = G4KaonZeroShort::KaonZeroShort(); 234 } 235 } 236 237 // KL + nucleus 238 else if (pdg == 130) { 239 G4double prob = (G4double)Z/(G4double)A; 240 if (G4UniformRand() >= prob) { 241 pd = G4KaonMinus::KaonMinus(); 242 } 243 else { 244 pd = G4KaonPlus::KaonPlus(); 245 } 246 } 247 248 return pd; 249 } 250 251 G4double 252 G4ChargeExchangeXS::ComputeDeuteronFraction(co 253 { 254 for (auto const & elm : *mat->GetElementVect 255 if (1 == elm->GetZasInt()) { 256 G4double ab = 0.0; 257 const G4int nIso = (G4int)elm->GetNumber 258 const G4double* abu = elm->GetRelativeAb 259 for (G4int j = 0; j < nIso; ++j) { 260 auto const iso = elm->GetIsotope(j); 261 ab += (iso->GetN() - iso->GetZ())*abu[j]; 262 } 263 return ab; 264 } 265 } 266 return 0.0; 267 } 268