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