| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // G4MuonicAtomHelper class implementation
27 //
28 // Author: K.Lynch, 01.07.2016 - First implementation
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::ConstructMuonicAtom(const G4String& name, G4int encoding,
43 G4Ions const* baseion)
44 {
45 // what should static charge be? for G4Ions, it is Z ... should it
46 // be Z-1 here (since there will always be a muon attached), or Z?
47 const G4double charge = baseion->GetPDGCharge();
48
49 static const G4String pType("MuonicAtom"); // put in an include? in an enum?
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::GetParticleTable()->FindParticle("mu-"))->GetPDGMass()
61 + baseion->GetPDGMass() - GetKShellEnergy(G4double(Z)); // fixme check
62
63 auto decayTable = new G4DecayTable();
64 // clang-format off
65 auto muatom = new G4MuonicAtom(name, mass, 0.0, charge,
66 baseion->GetPDGiSpin(),
67 baseion->GetPDGiParity(),
68 baseion->GetPDGiConjugation(),
69 baseion->GetPDGiIsospin(),
70 baseion->GetPDGiIsospin3(),
71 baseion->GetPDGiGParity(),
72 pType,
73 baseion->GetLeptonNumber(),
74 baseion->GetBaryonNumber(),
75 encoding,
76 stable,
77 lifetime,
78 decayTable,
79 shortlived,
80 baseion->GetParticleSubType(),
81 baseion);
82 // clang-format on
83
84 muatom->SetPDGMagneticMoment(baseion->GetPDGMagneticMoment());
85
86 // by this time both G4MuonicAtom and baseion should exist
87
88 // if we choose DIO this will be the channel we'll use, so we put
89 // br=1. to force it in that case
90
91 decayTable->Insert(new G4PhaseSpaceDecayChannel(name, 1., 4, "e-", "anti_nu_e", "nu_mu",
92 baseion->GetParticleName()));
93
94 muatom->SetDIOLifeTime(1. / lambdad);
95 muatom->SetNCLifeTime(1. / lambdac);
96 return muatom;
97 }
98
99 G4double G4MuonicAtomHelper::GetMuonCaptureRate(G4int Z, G4int A)
100 {
101 // Initialize data
102
103 // Mu- capture data from
104 // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
105 // weighted average of the two most precise measurements
106
107 // Data for Hydrogen from Phys. Rev. Lett. 99(2007)032002
108 // Data for Helium from D.F. Measday Phys. Rep. 354(2001)243
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 initialized as we exit the
119 // loop once Z is above the stored value; cRErr are not used now but
120 // are included for completeness and future use
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) / sizeof(capRates[0]);
219 for (std::size_t j = 0; j < nCapRates; ++j) {
220 if (capRates[j].Z == Z && capRates[j].A == A) {
221 lambda = capRates[j].cRate / microsecond;
222 break;
223 }
224 // make sure the data is sorted for the next statement to work correctly
225 if (capRates[j].Z > Z) {
226 break;
227 }
228 }
229
230 if (lambda < 0.) {
231 // == Mu capture lifetime (Goulard and Primakoff PRC10(1974)2034.
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-> -9 suggested by user
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) / G4double(Z);
243 G4double r2 = 1.0 - xmu;
244 lambda = t1 * zeff2 * zeff2 * (r2 * r2) * (1.0 - (1.0 - xmu) * .75704)
245 * (a2ze * b0a + 1.0 - (a2ze - 1.0) * b0b
246 - G4double(2 * (A - Z) + std::abs(a2ze - 1.)) * b0c / G4double(A * 4));
247 }
248
249 return lambda;
250 }
251
252 G4double G4MuonicAtomHelper::GetMuonZeff(G4int Z)
253 {
254 // == Effective charges from
255 // "Total Nuclear Capture Rates for Negative Muons"
256 // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
257 // and if not present from
258 // Ford and Wills Nucl Phys 35(1962)295 or interpolated
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.61, 7.49, 8.32, 9.14,
265 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.13,20.66,21.12,21.61,
267 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.95,28.20,28.42,28.64,
269 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.76,31.90,32.05,32.19,
271 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.42,34.52,34.63,34.73,
273 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
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(G4int Z)
287 {
288 // Decay time on K-shell
289 // N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1.
290
291 // this is the "small Z" approximation formula (2.9)
292 // Lambda(bound)/Lambda(free) = 1-beta(Z*alpha)**2 with beta~=2.5
293 // we assume that Z is Zeff
294
295 // PDG 2012 muon lifetime value is 2.1969811(22) 10e-6s
296 // which when inverted gives 0.45517005 10e+6/s
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 initialized as we exit the
306 // loop once Z is above the stored value
307
308 constexpr decRate decRates[] = {{1, 0.4558514, 0.0000151}};
309
310 G4double lambda = -1.;
311
312 // size_t nDecRates = sizeof(decRates)/sizeof(decRates[0]);
313 // for (size_t j = 0; j < nDecRates; ++j) {
314 // if( decRates[j].Z == Z ) {
315 // lambda = decRates[j].dRate / microsecond;
316 // break;
317 // }
318 // // make sure the data is sorted for the next statement to work
319 // if (decRates[j].Z > Z) {break;}
320 // }
321
322 // we'll use the above code once we have more data
323 // since we only have one value we just assign it
324 if (Z == 1) {
325 lambda = decRates[0].dRate / microsecond;
326 }
327
328 if (lambda < 0.) {
329 constexpr G4double freeMuonDecayRate = 0.45517005 / microsecond;
330 lambda = 1.0;
331 G4double x = GetMuonZeff(Z) * fine_structure_const;
332 lambda -= 2.5 * x * x;
333 lambda *= freeMuonDecayRate;
334 }
335
336 return lambda;
337 }
338
339 // From G4MuMinusCaptureCascade
340 G4double G4MuonicAtomHelper::GetKShellEnergy(G4double Z)
341 {
342 // Calculate the Energy of K Mesoatom Level for this Element using
343 // the Energy of Hydrogen Atom taken into account finite size of the
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., 21., 24.,
349 26., 29., 32., 38., 40., 41., 44., 49., 53., 55.,
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.326,
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, ListZK, ListKEnergy, Z);
362
363 return KEnergy;
364 }
365
366 // From G4MuMinusCaptureCascade
367 G4double G4MuonicAtomHelper::GetLinApprox(G4int N, const G4double* const X, const G4double* const Y,
368 G4double Xuser)
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]) * (Xuser - X[i - 1]) / (X[i] - X[i - 1]);
385 }
386 return Yuser;
387 }
388