| 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 // G4MuonicAtomHelper class implementation 26 // G4MuonicAtomHelper class implementation
27 // 27 //
28 // Author: K.Lynch, 01.07.2016 - First impleme 28 // Author: K.Lynch, 01.07.2016 - First implementation
29 // Revision: 29 // Revision:
30 // - 12.06.2017, K L Genser 30 // - 12.06.2017, K L Genser
31 // ------------------------------------------- 31 // --------------------------------------------------------------------
32 32
33 #include "G4MuonicAtomHelper.hh" 33 #include "G4MuonicAtomHelper.hh"
34 <<
35 #include "G4DecayTable.hh" 34 #include "G4DecayTable.hh"
36 #include "G4ParticleTable.hh" <<
37 #include "G4PhaseSpaceDecayChannel.hh" 35 #include "G4PhaseSpaceDecayChannel.hh"
38 #include "G4PhysicalConstants.hh" << 36 #include "G4ParticleTable.hh"
39 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
>> 38 #include "G4PhysicalConstants.hh"
40 #include "Randomize.hh" 39 #include "Randomize.hh"
41 40
42 G4MuonicAtom* G4MuonicAtomHelper::ConstructMuo << 41 G4MuonicAtom* G4MuonicAtomHelper::
43 << 42 ConstructMuonicAtom(const G4String& name, G4int encoding, G4Ions const* baseion)
44 { 43 {
45 // what should static charge be? for G4Ions 44 // what should static charge be? for G4Ions, it is Z ... should it
46 // be Z-1 here (since there will always be a 45 // be Z-1 here (since there will always be a muon attached), or Z?
47 const G4double charge = baseion->GetPDGCharg 46 const G4double charge = baseion->GetPDGCharge();
48 47
49 static const G4String pType("MuonicAtom"); << 48 static const G4String pType("MuonicAtom"); // put in an include? in an enum?
50 49
51 G4bool stable = false; << 50 G4bool stable = false;
52 // Get the lifetime 51 // Get the lifetime
53 G4int Z = baseion->GetAtomicNumber(); 52 G4int Z = baseion->GetAtomicNumber();
54 G4int A = baseion->GetAtomicMass(); 53 G4int A = baseion->GetAtomicMass();
55 G4double lambdac = GetMuonCaptureRate(Z, A); 54 G4double lambdac = GetMuonCaptureRate(Z, A);
56 G4double lambdad = GetMuonDecayRate(Z); 55 G4double lambdad = GetMuonDecayRate(Z);
57 G4double lifetime = 1. / (lambdac + lambdad) << 56 G4double lifetime = 1./(lambdac+lambdad);
58 G4bool shortlived = false; << 57 G4bool shortlived = false;
59 58
60 const G4double mass = (G4ParticleTable::GetP << 59 const G4double mass =
61 + baseion->GetPDGMass( << 60 (G4ParticleTable::GetParticleTable()->FindParticle("mu-"))->GetPDGMass() +
>> 61 baseion->GetPDGMass() - GetKShellEnergy(G4double(Z)); //fixme check
62 62
63 auto decayTable = new G4DecayTable(); << 63 G4DecayTable* decayTable = new G4DecayTable();
64 // clang-format off <<
65 auto muatom = new G4MuonicAtom(name, mass, 0 64 auto muatom = new G4MuonicAtom(name, mass, 0.0, charge,
66 baseion->GetP 65 baseion->GetPDGiSpin(),
67 baseion->GetP 66 baseion->GetPDGiParity(),
68 baseion->GetP 67 baseion->GetPDGiConjugation(),
69 baseion->GetP 68 baseion->GetPDGiIsospin(),
70 baseion->GetP 69 baseion->GetPDGiIsospin3(),
71 baseion->GetP 70 baseion->GetPDGiGParity(),
72 pType, 71 pType,
73 baseion->GetL 72 baseion->GetLeptonNumber(),
74 baseion->GetB 73 baseion->GetBaryonNumber(),
75 encoding, 74 encoding,
76 stable, 75 stable,
77 lifetime, 76 lifetime,
78 decayTable, 77 decayTable,
79 shortlived, 78 shortlived,
80 baseion->GetP 79 baseion->GetParticleSubType(),
81 baseion); 80 baseion);
82 // clang-format on <<
83 81
84 muatom->SetPDGMagneticMoment(baseion->GetPDG 82 muatom->SetPDGMagneticMoment(baseion->GetPDGMagneticMoment());
85 83
86 // by this time both G4MuonicAtom and baseio 84 // by this time both G4MuonicAtom and baseion should exist
87 85
88 // if we choose DIO this will be the channel 86 // if we choose DIO this will be the channel we'll use, so we put
89 // br=1. to force it in that case 87 // br=1. to force it in that case
90 88
91 decayTable->Insert(new G4PhaseSpaceDecayChan << 89 decayTable->Insert(new G4PhaseSpaceDecayChannel(name, 1., 4,
>> 90 "e-","anti_nu_e","nu_mu",
92 91 baseion->GetParticleName()));
93 92
94 muatom->SetDIOLifeTime(1. / lambdad); << 93 muatom->SetDIOLifeTime(1./lambdad);
95 muatom->SetNCLifeTime(1. / lambdac); << 94 muatom->SetNCLifeTime(1./lambdac);
96 return muatom; 95 return muatom;
>> 96
97 } 97 }
98 98
99 G4double G4MuonicAtomHelper::GetMuonCaptureRat 99 G4double G4MuonicAtomHelper::GetMuonCaptureRate(G4int Z, G4int A)
100 { 100 {
>> 101
101 // Initialize data 102 // Initialize data
102 103
103 // Mu- capture data from 104 // Mu- capture data from
104 // T. Suzuki, D. F. Measday, J.P. Roalsvig P 105 // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
105 // weighted average of the two most precise 106 // weighted average of the two most precise measurements
106 107
107 // Data for Hydrogen from Phys. Rev. Lett. 9 108 // Data for Hydrogen from Phys. Rev. Lett. 99(2007)032002
108 // Data for Helium from D.F. Measday Phys. R 109 // Data for Helium from D.F. Measday Phys. Rep. 354(2001)243
109 110
110 struct capRate << 111 struct capRate {
111 { << 112 G4int Z;
112 G4int Z; << 113 G4int A;
113 G4int A; << 114 G4double cRate;
114 G4double cRate; << 115 G4double cRErr;
115 G4double cRErr; <<
116 }; 116 };
117 117
118 // this struct has to be sorted by Z when in 118 // this struct has to be sorted by Z when initialized as we exit the
119 // loop once Z is above the stored value; cR 119 // loop once Z is above the stored value; cRErr are not used now but
120 // are included for completeness and future 120 // are included for completeness and future use
121 // clang-format off << 121
122 constexpr capRate capRates [] = { 122 constexpr capRate capRates [] = {
123 { 1, 1, 0.000725, 0.000017 }, 123 { 1, 1, 0.000725, 0.000017 },
124 { 2, 3, 0.002149, 0.00017 }, 124 { 2, 3, 0.002149, 0.00017 },
125 { 2, 4, 0.000356, 0.000026 }, 125 { 2, 4, 0.000356, 0.000026 },
126 { 3, 6, 0.004647, 0.00012 }, 126 { 3, 6, 0.004647, 0.00012 },
127 { 3, 7, 0.002229, 0.00012 }, 127 { 3, 7, 0.002229, 0.00012 },
128 { 4, 9, 0.006107, 0.00019 }, 128 { 4, 9, 0.006107, 0.00019 },
129 { 5, 10, 0.02757 , 0.00063 }, 129 { 5, 10, 0.02757 , 0.00063 },
130 { 5, 11, 0.02188 , 0.00064 }, 130 { 5, 11, 0.02188 , 0.00064 },
131 { 6, 12, 0.03807 , 0.00031 }, 131 { 6, 12, 0.03807 , 0.00031 },
132 { 6, 13, 0.03474 , 0.00034 }, 132 { 6, 13, 0.03474 , 0.00034 },
133 { 7, 14, 0.06885 , 0.00057 }, 133 { 7, 14, 0.06885 , 0.00057 },
134 { 8, 16, 0.10242 , 0.00059 }, 134 { 8, 16, 0.10242 , 0.00059 },
135 { 8, 18, 0.0880 , 0.0015 }, 135 { 8, 18, 0.0880 , 0.0015 },
136 { 9, 19, 0.22905 , 0.00099 }, 136 { 9, 19, 0.22905 , 0.00099 },
137 { 10, 20, 0.2288 , 0.0045 }, 137 { 10, 20, 0.2288 , 0.0045 },
138 { 11, 23, 0.3773 , 0.0014 }, 138 { 11, 23, 0.3773 , 0.0014 },
139 { 12, 24, 0.4823 , 0.0013 }, 139 { 12, 24, 0.4823 , 0.0013 },
140 { 13, 27, 0.6985 , 0.0012 }, 140 { 13, 27, 0.6985 , 0.0012 },
141 { 14, 28, 0.8656 , 0.0015 }, 141 { 14, 28, 0.8656 , 0.0015 },
142 { 15, 31, 1.1681 , 0.0026 }, 142 { 15, 31, 1.1681 , 0.0026 },
143 { 16, 32, 1.3510 , 0.0029 }, 143 { 16, 32, 1.3510 , 0.0029 },
144 { 17, 35, 1.800 , 0.050 }, 144 { 17, 35, 1.800 , 0.050 },
145 { 17, 37, 1.250 , 0.050 }, 145 { 17, 37, 1.250 , 0.050 },
146 { 18, 40, 1.2727 , 0.0650 }, 146 { 18, 40, 1.2727 , 0.0650 },
147 { 19, 39, 1.8492 , 0.0050 }, 147 { 19, 39, 1.8492 , 0.0050 },
148 { 20, 40, 2.5359 , 0.0070 }, 148 { 20, 40, 2.5359 , 0.0070 },
149 { 21, 45, 2.711 , 0.025 }, 149 { 21, 45, 2.711 , 0.025 },
150 { 22, 48, 2.5908 , 0.0115 }, 150 { 22, 48, 2.5908 , 0.0115 },
151 { 23, 51, 3.073 , 0.022 }, 151 { 23, 51, 3.073 , 0.022 },
152 { 24, 50, 3.825 , 0.050 }, 152 { 24, 50, 3.825 , 0.050 },
153 { 24, 52, 3.465 , 0.026 }, 153 { 24, 52, 3.465 , 0.026 },
154 { 24, 53, 3.297 , 0.045 }, 154 { 24, 53, 3.297 , 0.045 },
155 { 24, 54, 3.057 , 0.042 }, 155 { 24, 54, 3.057 , 0.042 },
156 { 25, 55, 3.900 , 0.030 }, 156 { 25, 55, 3.900 , 0.030 },
157 { 26, 56, 4.408 , 0.022 }, 157 { 26, 56, 4.408 , 0.022 },
158 { 27, 59, 4.945 , 0.025 }, 158 { 27, 59, 4.945 , 0.025 },
159 { 28, 58, 6.11 , 0.10 }, 159 { 28, 58, 6.11 , 0.10 },
160 { 28, 60, 5.56 , 0.10 }, 160 { 28, 60, 5.56 , 0.10 },
161 { 28, 62, 4.72 , 0.10 }, 161 { 28, 62, 4.72 , 0.10 },
162 { 29, 63, 5.691 , 0.030 }, 162 { 29, 63, 5.691 , 0.030 },
163 { 30, 66, 5.806 , 0.031 }, 163 { 30, 66, 5.806 , 0.031 },
164 { 31, 69, 5.700 , 0.060 }, 164 { 31, 69, 5.700 , 0.060 },
165 { 32, 72, 5.561 , 0.031 }, 165 { 32, 72, 5.561 , 0.031 },
166 { 33, 75, 6.094 , 0.037 }, 166 { 33, 75, 6.094 , 0.037 },
167 { 34, 80, 5.687 , 0.030 }, 167 { 34, 80, 5.687 , 0.030 },
168 { 35, 79, 7.223 , 0.28 }, 168 { 35, 79, 7.223 , 0.28 },
169 { 35, 81, 7.547 , 0.48 }, 169 { 35, 81, 7.547 , 0.48 },
170 { 37, 85, 6.89 , 0.14 }, 170 { 37, 85, 6.89 , 0.14 },
171 { 38, 88, 6.93 , 0.12 }, 171 { 38, 88, 6.93 , 0.12 },
172 { 39, 89, 7.89 , 0.11 }, 172 { 39, 89, 7.89 , 0.11 },
173 { 40, 91, 8.620 , 0.053 }, 173 { 40, 91, 8.620 , 0.053 },
174 { 41, 93, 10.38 , 0.11 }, 174 { 41, 93, 10.38 , 0.11 },
175 { 42, 96, 9.298 , 0.063 }, 175 { 42, 96, 9.298 , 0.063 },
176 { 45, 103, 10.010 , 0.045 }, 176 { 45, 103, 10.010 , 0.045 },
177 { 46, 106, 10.000 , 0.070 }, 177 { 46, 106, 10.000 , 0.070 },
178 { 47, 107, 10.869 , 0.095 }, 178 { 47, 107, 10.869 , 0.095 },
179 { 48, 112, 10.624 , 0.094 }, 179 { 48, 112, 10.624 , 0.094 },
180 { 49, 115, 11.38 , 0.11 }, 180 { 49, 115, 11.38 , 0.11 },
181 { 50, 119, 10.60 , 0.11 }, 181 { 50, 119, 10.60 , 0.11 },
182 { 51, 121, 10.40 , 0.12 }, 182 { 51, 121, 10.40 , 0.12 },
183 { 52, 128, 9.174 , 0.074 }, 183 { 52, 128, 9.174 , 0.074 },
184 { 53, 127, 11.276 , 0.098 }, 184 { 53, 127, 11.276 , 0.098 },
185 { 55, 133, 10.98 , 0.25 }, 185 { 55, 133, 10.98 , 0.25 },
186 { 56, 138, 10.112 , 0.085 }, 186 { 56, 138, 10.112 , 0.085 },
187 { 57, 139, 10.71 , 0.10 }, 187 { 57, 139, 10.71 , 0.10 },
188 { 58, 140, 11.501 , 0.087 }, 188 { 58, 140, 11.501 , 0.087 },
189 { 59, 141, 13.45 , 0.13 }, 189 { 59, 141, 13.45 , 0.13 },
190 { 60, 144, 12.35 , 0.13 }, 190 { 60, 144, 12.35 , 0.13 },
191 { 62, 150, 12.22 , 0.17 }, 191 { 62, 150, 12.22 , 0.17 },
192 { 64, 157, 12.00 , 0.13 }, 192 { 64, 157, 12.00 , 0.13 },
193 { 65, 159, 12.73 , 0.13 }, 193 { 65, 159, 12.73 , 0.13 },
194 { 66, 163, 12.29 , 0.18 }, 194 { 66, 163, 12.29 , 0.18 },
195 { 67, 165, 12.95 , 0.13 }, 195 { 67, 165, 12.95 , 0.13 },
196 { 68, 167, 13.04 , 0.27 }, 196 { 68, 167, 13.04 , 0.27 },
197 { 72, 178, 13.03 , 0.21 }, 197 { 72, 178, 13.03 , 0.21 },
198 { 73, 181, 12.86 , 0.13 }, 198 { 73, 181, 12.86 , 0.13 },
199 { 74, 184, 12.76 , 0.16 }, 199 { 74, 184, 12.76 , 0.16 },
200 { 79, 197, 13.35 , 0.10 }, 200 { 79, 197, 13.35 , 0.10 },
201 { 80, 201, 12.74 , 0.18 }, 201 { 80, 201, 12.74 , 0.18 },
202 { 81, 205, 13.85 , 0.17 }, 202 { 81, 205, 13.85 , 0.17 },
203 { 82, 207, 13.295 , 0.071 }, 203 { 82, 207, 13.295 , 0.071 },
204 { 83, 209, 13.238 , 0.065 }, 204 { 83, 209, 13.238 , 0.065 },
205 { 90, 232, 12.555 , 0.049 }, 205 { 90, 232, 12.555 , 0.049 },
206 { 92, 238, 12.592 , 0.035 }, 206 { 92, 238, 12.592 , 0.035 },
207 { 92, 233, 14.27 , 0.15 }, 207 { 92, 233, 14.27 , 0.15 },
208 { 92, 235, 13.470 , 0.085 }, 208 { 92, 235, 13.470 , 0.085 },
209 { 92, 236, 13.90 , 0.40 }, 209 { 92, 236, 13.90 , 0.40 },
210 { 93, 237, 13.58 , 0.18 }, 210 { 93, 237, 13.58 , 0.18 },
211 { 94, 239, 13.90 , 0.20 }, 211 { 94, 239, 13.90 , 0.20 },
212 { 94, 242, 12.86 , 0.19 } 212 { 94, 242, 12.86 , 0.19 }
213 }; 213 };
214 // clang-format on <<
215 214
216 G4double lambda = -1.; 215 G4double lambda = -1.;
217 216
218 std::size_t nCapRates = sizeof(capRates) / s << 217 std::size_t nCapRates = sizeof(capRates)/sizeof(capRates[0]);
219 for (std::size_t j = 0; j < nCapRates; ++j) << 218 for (std::size_t j = 0; j < nCapRates; ++j)
220 if (capRates[j].Z == Z && capRates[j].A == << 219 {
>> 220 if( capRates[j].Z == Z && capRates[j].A == A )
>> 221 {
221 lambda = capRates[j].cRate / microsecond 222 lambda = capRates[j].cRate / microsecond;
222 break; 223 break;
223 } 224 }
224 // make sure the data is sorted for the ne 225 // make sure the data is sorted for the next statement to work correctly
225 if (capRates[j].Z > Z) { << 226 if (capRates[j].Z > Z) {break;}
226 break; <<
227 } <<
228 } 227 }
229 228
230 if (lambda < 0.) { << 229 if (lambda < 0.)
>> 230 {
231 // == Mu capture lifetime (Goulard and Pr 231 // == Mu capture lifetime (Goulard and Primakoff PRC10(1974)2034.
232 232
233 constexpr G4double b0a = -0.03; 233 constexpr G4double b0a = -0.03;
234 constexpr G4double b0b = -0.25; 234 constexpr G4double b0b = -0.25;
235 constexpr G4double b0c = 3.24; << 235 constexpr G4double b0c = 3.24;
236 constexpr G4double t1 = 875.e-9; // -10-> << 236 constexpr G4double t1 = 875.e-9; // -10-> -9 suggested by user
237 G4double r1 = GetMuonZeff(Z); 237 G4double r1 = GetMuonZeff(Z);
238 G4double zeff2 = r1 * r1; 238 G4double zeff2 = r1 * r1;
239 239
240 // ^-4 -> ^-5 suggested by user 240 // ^-4 -> ^-5 suggested by user
241 G4double xmu = zeff2 * 2.663e-5; 241 G4double xmu = zeff2 * 2.663e-5;
242 G4double a2ze = 0.5 * G4double(A) / G4doub << 242 G4double a2ze = 0.5 *G4double(A) / G4double(Z);
243 G4double r2 = 1.0 - xmu; 243 G4double r2 = 1.0 - xmu;
244 lambda = t1 * zeff2 * zeff2 * (r2 * r2) * << 244 lambda = t1 * zeff2 * zeff2 * (r2 * r2) * (1.0 - (1.0 - xmu) * .75704) *
245 * (a2ze * b0a + 1.0 - (a2ze - 1.0 << 245 (a2ze * b0a + 1.0 - (a2ze - 1.0) * b0b -
246 - G4double(2 * (A - Z) + std:: << 246 G4double(2 * (A - Z) + std::abs(a2ze - 1.) ) * b0c / G4double(A * 4) );
>> 247
247 } 248 }
248 249
249 return lambda; 250 return lambda;
>> 251
250 } 252 }
251 253
252 G4double G4MuonicAtomHelper::GetMuonZeff(G4int 254 G4double G4MuonicAtomHelper::GetMuonZeff(G4int Z)
253 { 255 {
254 // == Effective charges from 256 // == Effective charges from
255 // "Total Nuclear Capture Rates for Negative 257 // "Total Nuclear Capture Rates for Negative Muons"
256 // T. Suzuki, D. F. Measday, J.P. Roalsvig P 258 // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
257 // and if not present from 259 // and if not present from
258 // Ford and Wills Nucl Phys 35(1962)295 or i 260 // Ford and Wills Nucl Phys 35(1962)295 or interpolated
259 261
260 // clang-format off <<
261 constexpr size_t maxZ = 100; 262 constexpr size_t maxZ = 100;
262 constexpr G4double zeff[maxZ+1] = 263 constexpr G4double zeff[maxZ+1] =
263 { 0., 264 { 0.,
264 1.00, 1.98, 2.94, 3.89, 4.81, 5.72, 6.6 265 1.00, 1.98, 2.94, 3.89, 4.81, 5.72, 6.61, 7.49, 8.32, 9.14,
265 9.95,10.69,11.48,12.22,12.90,13.64,14.2 266 9.95,10.69,11.48,12.22,12.90,13.64,14.24,14.89,15.53,16.15,
266 16.77,17.38,18.04,18.49,19.06,19.59,20. 267 16.77,17.38,18.04,18.49,19.06,19.59,20.13,20.66,21.12,21.61,
267 22.02,22.43,22.84,23.24,23.65,24.06,24. 268 22.02,22.43,22.84,23.24,23.65,24.06,24.47,24.85,25.23,25.61,
268 25.99,26.37,26.69,27.00,27.32,27.63,27. 269 25.99,26.37,26.69,27.00,27.32,27.63,27.95,28.20,28.42,28.64,
269 28.79,29.03,29.27,29.51,29.75,29.99,30. 270 28.79,29.03,29.27,29.51,29.75,29.99,30.22,30.36,30.53,30.69,
270 30.85,31.01,31.18,31.34,31.48,31.62,31. 271 30.85,31.01,31.18,31.34,31.48,31.62,31.76,31.90,32.05,32.19,
271 32.33,32.47,32.61,32.76,32.94,33.11,33. 272 32.33,32.47,32.61,32.76,32.94,33.11,33.29,33.46,33.64,33.81,
272 34.21,34.18,34.00,34.10,34.21,34.31,34. 273 34.21,34.18,34.00,34.10,34.21,34.31,34.42,34.52,34.63,34.73,
273 34.84,34.94,35.05,35.16,35.25,35.36,35. 274 34.84,34.94,35.05,35.16,35.25,35.36,35.46,35.57,35.67,35.78 };
274 // clang-format on <<
275 275
276 if (Z < 0) { << 276 if (Z<0) {Z=0;}
277 Z = 0; << 277 if (Z>G4int(maxZ)) {Z=maxZ;}
278 } <<
279 if (Z > G4int(maxZ)) { <<
280 Z = maxZ; <<
281 } <<
282 278
283 return zeff[Z]; 279 return zeff[Z];
284 } 280 }
285 281
286 G4double G4MuonicAtomHelper::GetMuonDecayRate( 282 G4double G4MuonicAtomHelper::GetMuonDecayRate(G4int Z)
287 { 283 {
288 // Decay time on K-shell 284 // Decay time on K-shell
289 // N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1. 285 // N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1.
290 286
291 // this is the "small Z" approximation formu 287 // this is the "small Z" approximation formula (2.9)
292 // Lambda(bound)/Lambda(free) = 1-beta(Z*alp 288 // Lambda(bound)/Lambda(free) = 1-beta(Z*alpha)**2 with beta~=2.5
293 // we assume that Z is Zeff 289 // we assume that Z is Zeff
294 290
295 // PDG 2012 muon lifetime value is 2.1969811 291 // PDG 2012 muon lifetime value is 2.1969811(22) 10e-6s
296 // which when inverted gives 0.4551700 292 // which when inverted gives 0.45517005 10e+6/s
297 293
298 struct decRate << 294 struct decRate {
299 { << 295 G4int Z;
300 G4int Z; << 296 G4double dRate;
301 G4double dRate; << 297 G4double dRErr;
302 G4double dRErr; <<
303 }; 298 };
304 299
305 // this struct has to be sorted by Z when in 300 // this struct has to be sorted by Z when initialized as we exit the
306 // loop once Z is above the stored value 301 // loop once Z is above the stored value
307 302
308 constexpr decRate decRates[] = {{1, 0.455851 << 303 constexpr decRate decRates [] = {
>> 304 { 1, 0.4558514, 0.0000151 }
>> 305 };
309 306
310 G4double lambda = -1.; 307 G4double lambda = -1.;
311 308
312 // size_t nDecRates = sizeof(decRates)/sizeo 309 // size_t nDecRates = sizeof(decRates)/sizeof(decRates[0]);
313 // for (size_t j = 0; j < nDecRates; ++j) { 310 // for (size_t j = 0; j < nDecRates; ++j) {
314 // if( decRates[j].Z == Z ) { 311 // if( decRates[j].Z == Z ) {
315 // lambda = decRates[j].dRate / microsec 312 // lambda = decRates[j].dRate / microsecond;
316 // break; 313 // break;
317 // } 314 // }
318 // // make sure the data is sorted for the 315 // // make sure the data is sorted for the next statement to work
319 // if (decRates[j].Z > Z) {break;} 316 // if (decRates[j].Z > Z) {break;}
320 // } 317 // }
321 318
322 // we'll use the above code once we have mor 319 // we'll use the above code once we have more data
323 // since we only have one value we just assi 320 // since we only have one value we just assign it
324 if (Z == 1) { << 321 if (Z == 1) { lambda = decRates[0].dRate/microsecond; }
325 lambda = decRates[0].dRate / microsecond; <<
326 } <<
327 322
328 if (lambda < 0.) { << 323 if (lambda < 0.)
329 constexpr G4double freeMuonDecayRate = 0.4 << 324 {
>> 325 constexpr G4double freeMuonDecayRate = 0.45517005 / microsecond;
330 lambda = 1.0; 326 lambda = 1.0;
331 G4double x = GetMuonZeff(Z) * fine_structu << 327 G4double x = GetMuonZeff(Z)*fine_structure_const;
332 lambda -= 2.5 * x * x; 328 lambda -= 2.5 * x * x;
333 lambda *= freeMuonDecayRate; 329 lambda *= freeMuonDecayRate;
334 } 330 }
335 331
336 return lambda; 332 return lambda;
337 } 333 }
338 334
339 // From G4MuMinusCaptureCascade 335 // From G4MuMinusCaptureCascade
340 G4double G4MuonicAtomHelper::GetKShellEnergy(G 336 G4double G4MuonicAtomHelper::GetKShellEnergy(G4double Z)
341 { 337 {
342 // Calculate the Energy of K Mesoatom Level 338 // Calculate the Energy of K Mesoatom Level for this Element using
343 // the Energy of Hydrogen Atom taken into ac 339 // the Energy of Hydrogen Atom taken into account finite size of the
344 // nucleus (V.Ivanchenko) 340 // nucleus (V.Ivanchenko)
345 // clang-format off <<
346 constexpr G4int ListK = 28; 341 constexpr G4int ListK = 28;
347 constexpr G4double ListZK[ListK] = { 342 constexpr G4double ListZK[ListK] = {
348 1., 2., 4., 6., 8., 11., 14., 17., 18 343 1., 2., 4., 6., 8., 11., 14., 17., 18., 21., 24.,
349 26., 29., 32., 38., 40., 41., 44., 49., 5 344 26., 29., 32., 38., 40., 41., 44., 49., 53., 55.,
350 60., 65., 70., 75., 81., 85., 92.}; 345 60., 65., 70., 75., 81., 85., 92.};
351 constexpr G4double ListKEnergy[ListK] = { 346 constexpr G4double ListKEnergy[ListK] = {
352 0.00275, 0.011, 0.043, 0.098, 0.173, 0.32 347 0.00275, 0.011, 0.043, 0.098, 0.173, 0.326,
353 0.524, 0.765, 0.853, 1.146, 1.472, 348 0.524, 0.765, 0.853, 1.146, 1.472,
354 1.708, 2.081, 2.475, 3.323, 3.627, 349 1.708, 2.081, 2.475, 3.323, 3.627,
355 3.779, 4.237, 5.016, 5.647, 5.966, 350 3.779, 4.237, 5.016, 5.647, 5.966,
356 6.793, 7.602, 8.421, 9.249, 10.222, 351 6.793, 7.602, 8.421, 9.249, 10.222,
357 10.923,11.984}; 352 10.923,11.984};
358 // clang-format on <<
359 353
360 // Energy with finite size corrections 354 // Energy with finite size corrections
361 G4double KEnergy = GetLinApprox(ListK, ListZ << 355 G4double KEnergy = GetLinApprox(ListK,ListZK,ListKEnergy,Z);
362 356
363 return KEnergy; 357 return KEnergy;
364 } 358 }
365 359
366 // From G4MuMinusCaptureCascade 360 // From G4MuMinusCaptureCascade
367 G4double G4MuonicAtomHelper::GetLinApprox(G4in << 361 G4double G4MuonicAtomHelper::GetLinApprox(G4int N,
>> 362 const G4double* const X,
>> 363 const G4double* const Y,
368 G4do 364 G4double Xuser)
369 { 365 {
370 G4double Yuser; 366 G4double Yuser;
371 if (Xuser <= X[0]) << 367 if(Xuser <= X[0]) Yuser = Y[0];
372 Yuser = Y[0]; << 368 else if(Xuser >= X[N-1]) Yuser = Y[N-1];
373 else if (Xuser >= X[N - 1]) << 369 else
374 Yuser = Y[N - 1]; << 370 {
375 else { <<
376 G4int i; 371 G4int i;
377 for (i = 1; i < N; ++i) { << 372 for (i=1; i<N; ++i)
378 if (Xuser <= X[i]) break; << 373 {
>> 374 if(Xuser <= X[i]) break;
379 } 375 }
380 376
381 if (Xuser == X[i]) << 377 if(Xuser == X[i]) Yuser = Y[i];
382 Yuser = Y[i]; << 378 else Yuser = Y[i-1] + (Y[i] - Y[i-1])*(Xuser - X[i-1])/(X[i] - X[i-1]);
383 else <<
384 Yuser = Y[i - 1] + (Y[i] - Y[i - 1]) * ( <<
385 } 379 }
386 return Yuser; 380 return Yuser;
387 } 381 }
388 382