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