Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: $ 26 // 27 // 27 // Author: D.H. Wright 28 // Author: D.H. Wright 28 // Date: 1 February 2011 29 // Date: 1 February 2011 29 // 30 // 30 // Modified: 31 // Modified: 31 // 32 // 32 // 19 Aug 2011, V.Ivanchenko move to new desig 33 // 19 Aug 2011, V.Ivanchenko move to new design and make x-section per element 33 34 34 // Description: use Kokoulin's parameterized c 35 // Description: use Kokoulin's parameterized calculation of virtual 35 // photon production cross sectio 36 // photon production cross section and conversion to 36 // real photons. 37 // real photons. 37 38 38 #include "G4KokoulinMuonNuclearXS.hh" 39 #include "G4KokoulinMuonNuclearXS.hh" 39 40 40 #include "G4PhysicalConstants.hh" 41 #include "G4PhysicalConstants.hh" 41 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh" >> 43 #include "G4PhysicsTable.hh" 42 #include "G4PhysicsLogVector.hh" 44 #include "G4PhysicsLogVector.hh" 43 #include "G4PhysicsVector.hh" << 44 #include "G4MuonMinus.hh" << 45 #include "G4MuonPlus.hh" << 46 #include "G4NucleiProperties.hh" << 47 #include "G4NistManager.hh" << 48 #include "G4Log.hh" << 49 #include "G4Exp.hh" << 50 45 51 // factory << 52 #include "G4CrossSectionFactory.hh" << 53 // << 54 G4_DECLARE_XS_FACTORY(G4KokoulinMuonNuclearXS) << 55 << 56 G4PhysicsVector* G4KokoulinMuonNuclearXS::theC << 57 46 58 G4KokoulinMuonNuclearXS::G4KokoulinMuonNuclear 47 G4KokoulinMuonNuclearXS::G4KokoulinMuonNuclearXS() 59 :G4VCrossSectionDataSet(Default_Name()), << 48 :G4VCrossSectionDataSet("KokoulinMuonNuclearXS"), theCrossSectionTable(0), 60 LowestKineticEnergy(1*GeV), HighestKineticEn 49 LowestKineticEnergy(1*GeV), HighestKineticEnergy(1*PeV), 61 TotBin(60), CutFixed(0.2*GeV), isInitialized << 50 TotBin(60), CutFixed(0.2*GeV) 62 {} << 63 << 64 G4KokoulinMuonNuclearXS::~G4KokoulinMuonNuclea << 65 { 51 { 66 if (isMaster) { << 52 BuildCrossSectionTable(); 67 for(G4int i=0; i<MAXZMUN; ++i) { << 68 delete theCrossSection[i]; << 69 theCrossSection[i] = 0; << 70 } << 71 } << 72 } 53 } 73 54 74 << 55 G4KokoulinMuonNuclearXS::~G4KokoulinMuonNuclearXS() 75 void << 76 G4KokoulinMuonNuclearXS::CrossSectionDescripti << 77 { 56 { 78 outFile << "G4KokoulinMuonNuclearXS provid << 57 if (theCrossSectionTable) { 79 << "cross section for mu- and mu+ interact << 58 theCrossSectionTable->clearAndDestroy(); 80 << "R. Kokoulin's approximation of the Bor << 59 delete theCrossSectionTable; 81 << "differential cross section at high ene << 60 } 82 << "over the muon energy loss to get the t << 83 << "function of muon kinetic energy\n" ; << 84 } 61 } 85 62 86 << 87 G4bool 63 G4bool 88 G4KokoulinMuonNuclearXS::IsElementApplicable(c 64 G4KokoulinMuonNuclearXS::IsElementApplicable(const G4DynamicParticle*, 89 G4int, const G4Material*) 65 G4int, const G4Material*) 90 { 66 { 91 return true; 67 return true; 92 } 68 } 93 69 94 void << 95 G4KokoulinMuonNuclearXS::BuildPhysicsTable(con << 96 { << 97 if(!isInitialized) { << 98 isInitialized = true; << 99 for(G4int i=0; i<MAXZMUN; ++i) { << 100 if(theCrossSection[i]) { return; } << 101 } << 102 isMaster = true; << 103 } << 104 if(isMaster) { BuildCrossSectionTable(); } << 105 } << 106 70 107 void G4KokoulinMuonNuclearXS::BuildCrossSectio 71 void G4KokoulinMuonNuclearXS::BuildCrossSectionTable() 108 { 72 { 109 G4double energy, A, Value; << 73 G4double lowEdgeEnergy; 110 G4int Z; << 74 G4PhysicsLogVector* ptrVector; 111 75 112 std::size_t nEl = G4Element::GetNumberOfElem << 76 G4int nEl = G4Element::GetNumberOfElements(); 113 const G4ElementTable* theElementTable = G4El 77 const G4ElementTable* theElementTable = G4Element::GetElementTable(); 114 G4NistManager* nistManager = G4NistManager:: << 115 78 116 for (std::size_t j = 0; j < nEl; ++j) { << 79 theCrossSectionTable = new G4PhysicsTable(nEl); 117 Z = G4lrint((*theElementTable)[j]->GetZ()) << 118 80 119 //AR-24Apr2018 Switch to treat transuranic << 81 G4double AtomicNumber; 120 const G4bool isHeavyElementAllowed = true; << 82 G4double AtomicWeight; 121 << 83 G4double Value; 122 A = nistManager->GetAtomicMassAmu(Z); << 84 123 if(Z < MAXZMUN && !theCrossSection[Z]) { << 85 for (G4int j = 0; j < nEl; j++) { 124 theCrossSection[Z] = new G4PhysicsLogVec << 86 ptrVector = new G4PhysicsLogVector(LowestKineticEnergy, 125 HighestKineticEnergy, << 87 HighestKineticEnergy, TotBin); 126 TotBin); << 88 AtomicNumber = (*theElementTable)[j]->GetZ(); 127 for (G4int i = 0; i <= TotBin; ++i) { << 89 AtomicWeight = (*theElementTable)[j]->GetA(); 128 energy = theCrossSection[Z]->Energy(i); << 90 129 Value = ComputeMicroscopicCrossSection(energ << 91 for (G4int i = 0; i <= TotBin; ++i) { 130 theCrossSection[Z]->PutValue(i,Value); << 92 lowEdgeEnergy = ptrVector->Energy(i); 131 } << 93 Value = ComputeMicroscopicCrossSection(lowEdgeEnergy, >> 94 AtomicNumber, AtomicWeight); >> 95 ptrVector->PutValue(i,Value); 132 } 96 } >> 97 >> 98 theCrossSectionTable->insertAt(j, ptrVector); >> 99 } >> 100 >> 101 // Build (Z,element) look-up table (for use in GetZandACrossSection) >> 102 for (G4int j = 0; j < nEl; j++) { >> 103 G4int ZZ = G4int((*theElementTable)[j]->GetZ()); >> 104 zelMap[ZZ] = j; 133 } 105 } 134 } 106 } 135 107 136 G4double G4KokoulinMuonNuclearXS:: 108 G4double G4KokoulinMuonNuclearXS:: 137 ComputeMicroscopicCrossSection(G4double Kineti << 109 ComputeMicroscopicCrossSection(G4double KineticEnergy, >> 110 G4double AtomicNumber, G4double AtomicWeight) 138 { 111 { 139 // Calculate cross section (differential in 112 // Calculate cross section (differential in muon incident kinetic energy) by 140 // integrating the double differential cross 113 // integrating the double differential cross section over the energy loss 141 114 142 static const G4double xgi[] = << 115 const G4double xgi[] = {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801}; 143 {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628 << 116 const G4double wgi[] = {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506}; 144 static const G4double wgi[] = << 117 const G4double ak1 = 6.9; 145 {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569 << 118 const G4double ak2 = 1.0; 146 static const G4double ak1 = 6.9; << 147 static const G4double ak2 = 1.0; << 148 119 149 G4double Mass = G4MuonMinus::MuonMinus()->Ge 120 G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass(); 150 121 151 G4double CrossSection = 0.0; 122 G4double CrossSection = 0.0; >> 123 if (AtomicNumber < 1.) return CrossSection; 152 if (KineticEnergy <= CutFixed) return CrossS 124 if (KineticEnergy <= CutFixed) return CrossSection; 153 125 154 G4double epmin = CutFixed; 126 G4double epmin = CutFixed; 155 G4double epmax = KineticEnergy + Mass - 0.5* 127 G4double epmax = KineticEnergy + Mass - 0.5*proton_mass_c2; 156 if (epmax <= epmin) return CrossSection; // 128 if (epmax <= epmin) return CrossSection; // NaN bug correction 157 129 158 G4double aaa = G4Log(epmin); << 130 G4double aaa = std::log(epmin) ; 159 G4double bbb = G4Log(epmax); << 131 G4double bbb = std::log(epmax) ; 160 G4int kkk = std::max(1,G4int((bbb-aaa)/ak1 + << 132 G4int kkk = G4int((bbb-aaa)/ak1 +ak2) ; 161 G4double hhh = (bbb-aaa)/kkk ; 133 G4double hhh = (bbb-aaa)/kkk ; 162 G4double epln; 134 G4double epln; 163 G4double ep; 135 G4double ep; 164 G4double x; 136 G4double x; 165 137 166 for (G4int l = 0; l < kkk; ++l) { << 138 for (G4int l = 0; l < kkk; l++) { 167 x = aaa + hhh*l; 139 x = aaa + hhh*l; 168 for (G4int ll = 0; ll < 8; ++ll) { << 140 for (G4int ll = 0; ll < 8; ll++) { 169 epln=x+xgi[ll]*hhh; 141 epln=x+xgi[ll]*hhh; 170 ep = G4Exp(epln); << 142 ep = std::exp(epln); 171 CrossSection += << 143 CrossSection += ep*wgi[ll]*ComputeDDMicroscopicCrossSection(KineticEnergy, 172 ep*wgi[ll]*ComputeDDMicroscopicCrossSection( << 144 AtomicNumber, AtomicWeight, ep); 173 } 145 } 174 } 146 } 175 147 176 CrossSection *= hhh ; 148 CrossSection *= hhh ; 177 if (CrossSection < 0.) { CrossSection = 0.; << 149 if (CrossSection < 0.) CrossSection = 0.; 178 return CrossSection; 150 return CrossSection; 179 } 151 } 180 152 181 G4double G4KokoulinMuonNuclearXS:: 153 G4double G4KokoulinMuonNuclearXS:: 182 ComputeDDMicroscopicCrossSection(G4double Kine << 154 ComputeDDMicroscopicCrossSection(G4double KineticEnergy, 183 G4double A, G << 155 G4double, G4double AtomicWeight, >> 156 G4double epsilon) 184 { 157 { 185 // Calculate the double-differential microsc 158 // Calculate the double-differential microscopic cross section (in muon 186 // incident kinetic energy and energy loss) 159 // incident kinetic energy and energy loss) using the cross section formula 187 // of R.P. Kokoulin (18/01/98) 160 // of R.P. Kokoulin (18/01/98) 188 161 189 static const G4double alam2 = 0.400*GeV*GeV; << 162 const G4double alam2 = 0.400*GeV*GeV; 190 static const G4double alam = 0.632456*GeV; << 163 const G4double alam = 0.632456*GeV; 191 static const G4double coeffn = fine_structur << 164 const G4double coeffn = fine_structure_const/pi; 192 165 193 G4double ParticleMass = G4MuonMinus::MuonMin 166 G4double ParticleMass = G4MuonMinus::MuonMinus()->GetPDGMass(); 194 G4double TotalEnergy = KineticEnergy + Parti 167 G4double TotalEnergy = KineticEnergy + ParticleMass; 195 168 196 G4double DCrossSection = 0.; 169 G4double DCrossSection = 0.; 197 170 198 if ((epsilon >= TotalEnergy - 0.5*proton_mas 171 if ((epsilon >= TotalEnergy - 0.5*proton_mass_c2) || 199 (epsilon <= CutFixed) ) { return DCrossS << 172 (epsilon <= CutFixed) ) return DCrossSection; 200 173 201 G4double ep = epsilon/GeV; 174 G4double ep = epsilon/GeV; 202 G4double aeff = 0.22*A+0.78*G4Exp(0.89*G4Log << 175 G4double a = AtomicWeight/(g/mole); 203 G4double sigph = (49.2+11.1*G4Log(ep)+151.8/ << 176 G4double aeff = 0.22*a+0.78*std::exp(0.89*std::log(a)); //shadowing >> 177 G4double sigph = (49.2+11.1*std::log(ep)+151.8/std::sqrt(ep))*microbarn; 204 178 205 G4double v = epsilon/TotalEnergy; 179 G4double v = epsilon/TotalEnergy; 206 G4double v1 = 1.-v; 180 G4double v1 = 1.-v; 207 G4double v2 = v*v; 181 G4double v2 = v*v; 208 G4double mass2 = ParticleMass*ParticleMass; 182 G4double mass2 = ParticleMass*ParticleMass; 209 183 210 G4double up = TotalEnergy*TotalEnergy*v1/mas 184 G4double up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1)); 211 G4double down = 1.+epsilon/alam*(1.+alam/(2. 185 G4double down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam); 212 186 213 DCrossSection = coeffn*aeff*sigph/epsilon* 187 DCrossSection = coeffn*aeff*sigph/epsilon* 214 (-v1+(v1+0.5*v2*(1.+2.*mass2 << 188 (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*std::log(up/down)); 215 189 216 if (DCrossSection < 0.) { DCrossSection = 0. << 190 if (DCrossSection < 0.) DCrossSection = 0.; 217 return DCrossSection; 191 return DCrossSection; 218 } 192 } 219 193 220 G4double G4KokoulinMuonNuclearXS:: 194 G4double G4KokoulinMuonNuclearXS:: 221 GetElementCrossSection(const G4DynamicParticle 195 GetElementCrossSection(const G4DynamicParticle* aPart, 222 G4int ZZ, const G4Material*) 196 G4int ZZ, const G4Material*) 223 { 197 { 224 //AR-24Apr2018 Switch to treat transuranic e << 198 G4int j = zelMap[ZZ]; 225 G4int Z = std::min(ZZ, 92); << 199 return (*theCrossSectionTable)[j]->Value(aPart->GetKineticEnergy()); 226 return theCrossSection[Z]->LogVectorValue(aP << 227 aPart->GetLogKineticEnergy()); << 228 } 200 } 229 201 230 202