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