| 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"
42 #include "G4PhysicsLogVector.hh" 43 #include "G4PhysicsLogVector.hh"
43 #include "G4PhysicsVector.hh" 44 #include "G4PhysicsVector.hh"
44 #include "G4MuonMinus.hh" 45 #include "G4MuonMinus.hh"
45 #include "G4MuonPlus.hh" 46 #include "G4MuonPlus.hh"
46 #include "G4NucleiProperties.hh" 47 #include "G4NucleiProperties.hh"
47 #include "G4NistManager.hh" 48 #include "G4NistManager.hh"
48 #include "G4Log.hh" 49 #include "G4Log.hh"
49 #include "G4Exp.hh" 50 #include "G4Exp.hh"
50 51
51 // factory 52 // factory
52 #include "G4CrossSectionFactory.hh" 53 #include "G4CrossSectionFactory.hh"
53 // 54 //
54 G4_DECLARE_XS_FACTORY(G4KokoulinMuonNuclearXS) 55 G4_DECLARE_XS_FACTORY(G4KokoulinMuonNuclearXS);
55 56
56 G4PhysicsVector* G4KokoulinMuonNuclearXS::theC 57 G4PhysicsVector* G4KokoulinMuonNuclearXS::theCrossSection[] = {0};
57 58
58 G4KokoulinMuonNuclearXS::G4KokoulinMuonNuclear 59 G4KokoulinMuonNuclearXS::G4KokoulinMuonNuclearXS()
59 :G4VCrossSectionDataSet(Default_Name()), 60 :G4VCrossSectionDataSet(Default_Name()),
60 LowestKineticEnergy(1*GeV), HighestKineticEn 61 LowestKineticEnergy(1*GeV), HighestKineticEnergy(1*PeV),
61 TotBin(60), CutFixed(0.2*GeV), isInitialized 62 TotBin(60), CutFixed(0.2*GeV), isInitialized(false), isMaster(false)
62 {} 63 {}
63 64
64 G4KokoulinMuonNuclearXS::~G4KokoulinMuonNuclea 65 G4KokoulinMuonNuclearXS::~G4KokoulinMuonNuclearXS()
65 { 66 {
66 if (isMaster) { 67 if (isMaster) {
67 for(G4int i=0; i<MAXZMUN; ++i) { 68 for(G4int i=0; i<MAXZMUN; ++i) {
68 delete theCrossSection[i]; 69 delete theCrossSection[i];
69 theCrossSection[i] = 0; 70 theCrossSection[i] = 0;
70 } 71 }
71 } 72 }
72 } 73 }
73 74
74 75
75 void 76 void
76 G4KokoulinMuonNuclearXS::CrossSectionDescripti 77 G4KokoulinMuonNuclearXS::CrossSectionDescription(std::ostream& outFile) const
77 { 78 {
78 outFile << "G4KokoulinMuonNuclearXS provid 79 outFile << "G4KokoulinMuonNuclearXS provides the total inelastic\n"
79 << "cross section for mu- and mu+ interact 80 << "cross section for mu- and mu+ interactions with nuclei.\n"
80 << "R. Kokoulin's approximation of the Bor 81 << "R. Kokoulin's approximation of the Borog and Petrukhin double\n"
81 << "differential cross section at high ene 82 << "differential cross section at high energy and low Q**2 is integrated\n"
82 << "over the muon energy loss to get the t 83 << "over the muon energy loss to get the total cross section as a\n"
83 << "function of muon kinetic energy\n" ; 84 << "function of muon kinetic energy\n" ;
84 } 85 }
85 86
86 87
87 G4bool 88 G4bool
88 G4KokoulinMuonNuclearXS::IsElementApplicable(c 89 G4KokoulinMuonNuclearXS::IsElementApplicable(const G4DynamicParticle*,
89 G4int, const G4Material*) 90 G4int, const G4Material*)
90 { 91 {
91 return true; 92 return true;
92 } 93 }
93 94
94 void 95 void
95 G4KokoulinMuonNuclearXS::BuildPhysicsTable(con 96 G4KokoulinMuonNuclearXS::BuildPhysicsTable(const G4ParticleDefinition&)
96 { 97 {
97 if(!isInitialized) { 98 if(!isInitialized) {
98 isInitialized = true; 99 isInitialized = true;
99 for(G4int i=0; i<MAXZMUN; ++i) { 100 for(G4int i=0; i<MAXZMUN; ++i) {
100 if(theCrossSection[i]) { return; } 101 if(theCrossSection[i]) { return; }
101 } 102 }
102 isMaster = true; 103 isMaster = true;
103 } 104 }
104 if(isMaster) { BuildCrossSectionTable(); } 105 if(isMaster) { BuildCrossSectionTable(); }
105 } 106 }
106 107
107 void G4KokoulinMuonNuclearXS::BuildCrossSectio 108 void G4KokoulinMuonNuclearXS::BuildCrossSectionTable()
108 { 109 {
109 G4double energy, A, Value; 110 G4double energy, A, Value;
110 G4int Z; 111 G4int Z;
111 112
112 std::size_t nEl = G4Element::GetNumberOfElem << 113 G4int nEl = G4Element::GetNumberOfElements();
113 const G4ElementTable* theElementTable = G4El 114 const G4ElementTable* theElementTable = G4Element::GetElementTable();
114 G4NistManager* nistManager = G4NistManager:: 115 G4NistManager* nistManager = G4NistManager::Instance();
115 116
116 for (std::size_t j = 0; j < nEl; ++j) { << 117 for (G4int j = 0; j < nEl; j++) {
117 Z = G4lrint((*theElementTable)[j]->GetZ()) 118 Z = G4lrint((*theElementTable)[j]->GetZ());
118 <<
119 //AR-24Apr2018 Switch to treat transuranic <<
120 const G4bool isHeavyElementAllowed = true; <<
121 <<
122 A = nistManager->GetAtomicMassAmu(Z); 119 A = nistManager->GetAtomicMassAmu(Z);
123 if(Z < MAXZMUN && !theCrossSection[Z]) { 120 if(Z < MAXZMUN && !theCrossSection[Z]) {
124 theCrossSection[Z] = new G4PhysicsLogVec 121 theCrossSection[Z] = new G4PhysicsLogVector(LowestKineticEnergy,
125 HighestKineticEnergy, 122 HighestKineticEnergy,
126 TotBin); 123 TotBin);
127 for (G4int i = 0; i <= TotBin; ++i) { 124 for (G4int i = 0; i <= TotBin; ++i) {
128 energy = theCrossSection[Z]->Energy(i); 125 energy = theCrossSection[Z]->Energy(i);
129 Value = ComputeMicroscopicCrossSection(energ 126 Value = ComputeMicroscopicCrossSection(energy, A);
130 theCrossSection[Z]->PutValue(i,Value); 127 theCrossSection[Z]->PutValue(i,Value);
131 } 128 }
132 } 129 }
133 } 130 }
134 } 131 }
135 132
136 G4double G4KokoulinMuonNuclearXS:: 133 G4double G4KokoulinMuonNuclearXS::
137 ComputeMicroscopicCrossSection(G4double Kineti 134 ComputeMicroscopicCrossSection(G4double KineticEnergy, G4double A)
138 { 135 {
139 // Calculate cross section (differential in 136 // Calculate cross section (differential in muon incident kinetic energy) by
140 // integrating the double differential cross 137 // integrating the double differential cross section over the energy loss
141 138
142 static const G4double xgi[] = 139 static const G4double xgi[] =
143 {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628 140 {0.0199,0.1017,0.2372,0.4083,0.5917,0.7628,0.8983,0.9801};
144 static const G4double wgi[] = 141 static const G4double wgi[] =
145 {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569 142 {0.0506,0.1112,0.1569,0.1813,0.1813,0.1569,0.1112,0.0506};
146 static const G4double ak1 = 6.9; 143 static const G4double ak1 = 6.9;
147 static const G4double ak2 = 1.0; 144 static const G4double ak2 = 1.0;
148 145
149 G4double Mass = G4MuonMinus::MuonMinus()->Ge 146 G4double Mass = G4MuonMinus::MuonMinus()->GetPDGMass();
150 147
151 G4double CrossSection = 0.0; 148 G4double CrossSection = 0.0;
152 if (KineticEnergy <= CutFixed) return CrossS 149 if (KineticEnergy <= CutFixed) return CrossSection;
153 150
154 G4double epmin = CutFixed; 151 G4double epmin = CutFixed;
155 G4double epmax = KineticEnergy + Mass - 0.5* 152 G4double epmax = KineticEnergy + Mass - 0.5*proton_mass_c2;
156 if (epmax <= epmin) return CrossSection; // 153 if (epmax <= epmin) return CrossSection; // NaN bug correction
157 154
158 G4double aaa = G4Log(epmin); 155 G4double aaa = G4Log(epmin);
159 G4double bbb = G4Log(epmax); 156 G4double bbb = G4Log(epmax);
160 G4int kkk = std::max(1,G4int((bbb-aaa)/ak1 + 157 G4int kkk = std::max(1,G4int((bbb-aaa)/ak1 +ak2));
161 G4double hhh = (bbb-aaa)/kkk ; 158 G4double hhh = (bbb-aaa)/kkk ;
162 G4double epln; 159 G4double epln;
163 G4double ep; 160 G4double ep;
164 G4double x; 161 G4double x;
165 162
166 for (G4int l = 0; l < kkk; ++l) { 163 for (G4int l = 0; l < kkk; ++l) {
167 x = aaa + hhh*l; 164 x = aaa + hhh*l;
168 for (G4int ll = 0; ll < 8; ++ll) { 165 for (G4int ll = 0; ll < 8; ++ll) {
169 epln=x+xgi[ll]*hhh; 166 epln=x+xgi[ll]*hhh;
170 ep = G4Exp(epln); 167 ep = G4Exp(epln);
171 CrossSection += 168 CrossSection +=
172 ep*wgi[ll]*ComputeDDMicroscopicCrossSection( 169 ep*wgi[ll]*ComputeDDMicroscopicCrossSection(KineticEnergy, 0, A, ep);
173 } 170 }
174 } 171 }
175 172
176 CrossSection *= hhh ; 173 CrossSection *= hhh ;
177 if (CrossSection < 0.) { CrossSection = 0.; 174 if (CrossSection < 0.) { CrossSection = 0.; }
178 return CrossSection; 175 return CrossSection;
179 } 176 }
180 177
181 G4double G4KokoulinMuonNuclearXS:: 178 G4double G4KokoulinMuonNuclearXS::
182 ComputeDDMicroscopicCrossSection(G4double Kine 179 ComputeDDMicroscopicCrossSection(G4double KineticEnergy, G4double,
183 G4double A, G 180 G4double A, G4double epsilon)
184 { 181 {
185 // Calculate the double-differential microsc 182 // Calculate the double-differential microscopic cross section (in muon
186 // incident kinetic energy and energy loss) 183 // incident kinetic energy and energy loss) using the cross section formula
187 // of R.P. Kokoulin (18/01/98) 184 // of R.P. Kokoulin (18/01/98)
188 185
189 static const G4double alam2 = 0.400*GeV*GeV; 186 static const G4double alam2 = 0.400*GeV*GeV;
190 static const G4double alam = 0.632456*GeV; 187 static const G4double alam = 0.632456*GeV;
191 static const G4double coeffn = fine_structur 188 static const G4double coeffn = fine_structure_const/pi;
192 189
193 G4double ParticleMass = G4MuonMinus::MuonMin 190 G4double ParticleMass = G4MuonMinus::MuonMinus()->GetPDGMass();
194 G4double TotalEnergy = KineticEnergy + Parti 191 G4double TotalEnergy = KineticEnergy + ParticleMass;
195 192
196 G4double DCrossSection = 0.; 193 G4double DCrossSection = 0.;
197 194
198 if ((epsilon >= TotalEnergy - 0.5*proton_mas 195 if ((epsilon >= TotalEnergy - 0.5*proton_mass_c2) ||
199 (epsilon <= CutFixed) ) { return DCrossS 196 (epsilon <= CutFixed) ) { return DCrossSection; }
200 197
201 G4double ep = epsilon/GeV; 198 G4double ep = epsilon/GeV;
202 G4double aeff = 0.22*A+0.78*G4Exp(0.89*G4Log 199 G4double aeff = 0.22*A+0.78*G4Exp(0.89*G4Log(A)); //shadowing
203 G4double sigph = (49.2+11.1*G4Log(ep)+151.8/ 200 G4double sigph = (49.2+11.1*G4Log(ep)+151.8/std::sqrt(ep))*microbarn;
204 201
205 G4double v = epsilon/TotalEnergy; 202 G4double v = epsilon/TotalEnergy;
206 G4double v1 = 1.-v; 203 G4double v1 = 1.-v;
207 G4double v2 = v*v; 204 G4double v2 = v*v;
208 G4double mass2 = ParticleMass*ParticleMass; 205 G4double mass2 = ParticleMass*ParticleMass;
209 206
210 G4double up = TotalEnergy*TotalEnergy*v1/mas 207 G4double up = TotalEnergy*TotalEnergy*v1/mass2*(1.+mass2*v2/(alam2*v1));
211 G4double down = 1.+epsilon/alam*(1.+alam/(2. 208 G4double down = 1.+epsilon/alam*(1.+alam/(2.*proton_mass_c2)+epsilon/alam);
212 209
213 DCrossSection = coeffn*aeff*sigph/epsilon* 210 DCrossSection = coeffn*aeff*sigph/epsilon*
214 (-v1+(v1+0.5*v2*(1.+2.*mass2 211 (-v1+(v1+0.5*v2*(1.+2.*mass2/alam2))*G4Log(up/down));
215 212
216 if (DCrossSection < 0.) { DCrossSection = 0. 213 if (DCrossSection < 0.) { DCrossSection = 0.; }
217 return DCrossSection; 214 return DCrossSection;
218 } 215 }
219 216
220 G4double G4KokoulinMuonNuclearXS:: 217 G4double G4KokoulinMuonNuclearXS::
221 GetElementCrossSection(const G4DynamicParticle 218 GetElementCrossSection(const G4DynamicParticle* aPart,
222 G4int ZZ, const G4Material*) << 219 G4int Z, const G4Material*)
223 { 220 {
224 //AR-24Apr2018 Switch to treat transuranic e << 221 return theCrossSection[Z]->Value(aPart->GetKineticEnergy());
225 G4int Z = std::min(ZZ, 92); <<
226 return theCrossSection[Z]->LogVectorValue(aP <<
227 aPart->GetLogKineticEnergy()); <<
228 } 222 }
229 223
230 224