ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MuonMinusBoundDecay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MuonMinusBoundDecay.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 // GEANT4 Class header file
30 //
31 // File name: G4MuonMinusBoundDecay
32 //
33 // Author: V.Ivanchenko (Vladimir.Ivantchenko@cern.ch)
34 //
35 // Creation date: 24 April 2012 on base of G4MuMinusCaptureAtRest
36 //
37 // Modified:
38 // 04/23/2013 K.Genser Fixed a constant in computation of lambda
39 // as suggested by J P Miller/Y Oksuzian;
40 // Optimized and corrected lambda calculation/lookup
41 // 04/30/2013 K.Genser Improved GetMuonCaptureRate extended data and lookup
42 // to take both Z & A into account
43 // Improved GetMuonDecayRate by using Zeff instead of Z
44 // Extracted Zeff into GetMuonZeff
45 // 10/08/2018 K.Genser Improved accuracy of the bound decay rate
46 //
47 //----------------------------------------------------------------------
48 
49 #include "G4MuonMinusBoundDecay.hh"
50 #include "Randomize.hh"
51 #include "G4RandomDirection.hh"
52 #include "G4PhysicalConstants.hh"
53 #include "G4SystemOfUnits.hh"
54 #include "G4ThreeVector.hh"
55 #include "G4NistManager.hh"
56 #include "G4NucleiProperties.hh"
57 #include "G4IonTable.hh"
58 #include "G4MuonMinus.hh"
59 #include "G4Electron.hh"
60 #include "G4NeutrinoMu.hh"
61 #include "G4AntiNeutrinoE.hh"
62 #include "G4Log.hh"
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
66 G4MuonMinusBoundDecay::G4MuonMinusBoundDecay()
67  : G4HadronicInteraction("muMinusBoundDeacy")
68 {
69  fMuMass = G4MuonMinus::MuonMinus()->GetPDGMass();
70 }
71 
72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73 
74 G4MuonMinusBoundDecay::~G4MuonMinusBoundDecay()
75 {}
76 
77 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
78 
79 G4HadFinalState*
80 G4MuonMinusBoundDecay::ApplyYourself(const G4HadProjectile& projectile,
81  G4Nucleus& targetNucleus)
82 {
83  result.Clear();
84  G4int Z = targetNucleus.GetZ_asInt();
85  G4int A = targetNucleus.GetA_asInt();
86 
87  // Decide on Decay or Capture, and doit.
88  G4double lambdac = GetMuonCaptureRate(Z, A);
89  G4double lambdad = GetMuonDecayRate(Z, A, fMuMass,
90  targetNucleus.AtomicMass(A,Z));
91  G4double lambda = lambdac + lambdad;
92 
93  // === sample capture time and change time of projectile
94  // === this is needed for the case when bound decay is not happen
95  // === but muon is capruted by the nucleus with some delay
96 
97  G4HadProjectile* p = const_cast<G4HadProjectile*>(&projectile);
98  G4double time = p->GetGlobalTime() - G4Log(G4UniformRand())/lambda;
99  p->SetGlobalTime(time);
100 
101  //G4cout << "lambda= " << lambda << " lambdac= " << lambdac
102  //<< " t= " << time << G4endl;
103 
104  // cascade
105  if( G4UniformRand()*lambda < lambdac) {
106  result.SetStatusChange(isAlive);
107 
108  } else {
109 
110  // Simulation on Decay of mu- on a K-shell of the muonic atom
111  result.SetStatusChange(stopAndKill);
112  G4double xmax = 1 + electron_mass_c2*electron_mass_c2/(fMuMass*fMuMass);
113  G4double xmin = 2.0*electron_mass_c2/fMuMass;
114  G4double KEnergy = projectile.GetBoundEnergy();
115 
116  /*
117  G4cout << "G4MuonMinusBoundDecay::ApplyYourself"
118  << " XMAX= " << xmax << " Ebound= " << KEnergy<< G4endl;
119  */
120  G4double pmu = std::sqrt(KEnergy*(KEnergy + 2.0*fMuMass));
121  G4double emu = KEnergy + fMuMass;
122  G4ThreeVector dir = G4RandomDirection();
123  G4LorentzVector MU(pmu*dir, emu);
124  G4ThreeVector bst = MU.boostVector();
125 
126  G4double Eelect, Pelect, x, ecm;
127  G4LorentzVector EL, NN;
128  // Calculate electron energy
129  // these do/while loops are safe
130  do {
131  do {
132  x = xmin + (xmax-xmin)*G4UniformRand();
133  } while (G4UniformRand() > (3.0 - 2.0*x)*x*x );
134  Eelect = x*fMuMass*0.5;
135  Pelect = 0.0;
136  if(Eelect > electron_mass_c2) {
137  Pelect = std::sqrt(Eelect*Eelect - electron_mass_c2*electron_mass_c2);
138  } else {
139  Pelect = 0.0;
140  Eelect = electron_mass_c2;
141  }
142  dir = G4RandomDirection();
143  EL = G4LorentzVector(Pelect*dir,Eelect);
144  EL.boost(bst);
145  Eelect = EL.e() - electron_mass_c2 - 2.0*KEnergy;
146  //
147  // Calculate rest frame parameters of 2 neutrinos
148  //
149  NN = MU - EL;
150  ecm = NN.mag2();
151  // Loop checking, 06-Aug-2015, Vladimir Ivanchenko
152  } while (Eelect < 0.0 || ecm < 0.0);
153 
154  //
155  // Create electron
156  //
157  G4DynamicParticle* dp = new G4DynamicParticle(G4Electron::Electron(),
158  EL.vect().unit(),
159  Eelect);
160 
161  AddNewParticle(dp, time);
162  //
163  // Create Neutrinos
164  //
165  ecm = 0.5*std::sqrt(ecm);
166  bst = NN.boostVector();
167  G4ThreeVector p1 = ecm * G4RandomDirection();
168  G4LorentzVector N1 = G4LorentzVector(p1,ecm);
169  N1.boost(bst);
170  dp = new G4DynamicParticle(G4AntiNeutrinoE::AntiNeutrinoE(), N1);
171  AddNewParticle(dp, time);
172  NN -= N1;
173  dp = new G4DynamicParticle(G4NeutrinoMu::NeutrinoMu(), NN);
174  AddNewParticle(dp, time);
175  }
176  return &result;
177 }
178 
179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
180 
181 G4double G4MuonMinusBoundDecay::GetMuonCaptureRate(G4int Z, G4int A)
182 {
183 
184  // Initialize data
185 
186  // Mu- capture data from
187  // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
188  // weighted average of the two most precise measurements
189 
190  // Data for Hydrogen from Phys. Rev. Lett. 99(2007)032002
191  // Data for Helium from D.F. Measday Phys. Rep. 354(2001)243
192 
193  struct capRate {
194  G4int Z;
195  G4int A;
196  G4double cRate;
197  G4double cRErr;
198  };
199 
200  // this struct has to be sorted by Z when initialized as we exit the
201  // loop once Z is above the stored value; cRErr are not used now but
202  // are included for completeness and future use
203 
204  static const capRate capRates [] = {
205  { 1, 1, 0.000725, 0.000017 },
206  { 2, 3, 0.002149, 0.00017 },
207  { 2, 4, 0.000356, 0.000026 },
208  { 3, 6, 0.004647, 0.00012 },
209  { 3, 7, 0.002229, 0.00012 },
210  { 4, 9, 0.006107, 0.00019 },
211  { 5, 10, 0.02757 , 0.00063 },
212  { 5, 11, 0.02188 , 0.00064 },
213  { 6, 12, 0.03807 , 0.00031 },
214  { 6, 13, 0.03474 , 0.00034 },
215  { 7, 14, 0.06885 , 0.00057 },
216  { 8, 16, 0.10242 , 0.00059 },
217  { 8, 18, 0.0880 , 0.0015 },
218  { 9, 19, 0.22905 , 0.00099 },
219  { 10, 20, 0.2288 , 0.0045 },
220  { 11, 23, 0.3773 , 0.0014 },
221  { 12, 24, 0.4823 , 0.0013 },
222  { 13, 27, 0.6985 , 0.0012 },
223  { 14, 28, 0.8656 , 0.0015 },
224  { 15, 31, 1.1681 , 0.0026 },
225  { 16, 32, 1.3510 , 0.0029 },
226  { 17, 35, 1.800 , 0.050 },
227  { 17, 37, 1.250 , 0.050 },
228  { 18, 40, 1.2727 , 0.0650 },
229  { 19, 39, 1.8492 , 0.0050 },
230  { 20, 40, 2.5359 , 0.0070 },
231  { 21, 45, 2.711 , 0.025 },
232  { 22, 48, 2.5908 , 0.0115 },
233  { 23, 51, 3.073 , 0.022 },
234  { 24, 50, 3.825 , 0.050 },
235  { 24, 52, 3.465 , 0.026 },
236  { 24, 53, 3.297 , 0.045 },
237  { 24, 54, 3.057 , 0.042 },
238  { 25, 55, 3.900 , 0.030 },
239  { 26, 56, 4.408 , 0.022 },
240  { 27, 59, 4.945 , 0.025 },
241  { 28, 58, 6.11 , 0.10 },
242  { 28, 60, 5.56 , 0.10 },
243  { 28, 62, 4.72 , 0.10 },
244  { 29, 63, 5.691 , 0.030 },
245  { 30, 66, 5.806 , 0.031 },
246  { 31, 69, 5.700 , 0.060 },
247  { 32, 72, 5.561 , 0.031 },
248  { 33, 75, 6.094 , 0.037 },
249  { 34, 80, 5.687 , 0.030 },
250  { 35, 79, 7.223 , 0.28 },
251  { 35, 81, 7.547 , 0.48 },
252  { 37, 85, 6.89 , 0.14 },
253  { 38, 88, 6.93 , 0.12 },
254  { 39, 89, 7.89 , 0.11 },
255  { 40, 91, 8.620 , 0.053 },
256  { 41, 93, 10.38 , 0.11 },
257  { 42, 96, 9.298 , 0.063 },
258  { 45, 103, 10.010 , 0.045 },
259  { 46, 106, 10.000 , 0.070 },
260  { 47, 107, 10.869 , 0.095 },
261  { 48, 112, 10.624 , 0.094 },
262  { 49, 115, 11.38 , 0.11 },
263  { 50, 119, 10.60 , 0.11 },
264  { 51, 121, 10.40 , 0.12 },
265  { 52, 128, 9.174 , 0.074 },
266  { 53, 127, 11.276 , 0.098 },
267  { 55, 133, 10.98 , 0.25 },
268  { 56, 138, 10.112 , 0.085 },
269  { 57, 139, 10.71 , 0.10 },
270  { 58, 140, 11.501 , 0.087 },
271  { 59, 141, 13.45 , 0.13 },
272  { 60, 144, 12.35 , 0.13 },
273  { 62, 150, 12.22 , 0.17 },
274  { 64, 157, 12.00 , 0.13 },
275  { 65, 159, 12.73 , 0.13 },
276  { 66, 163, 12.29 , 0.18 },
277  { 67, 165, 12.95 , 0.13 },
278  { 68, 167, 13.04 , 0.27 },
279  { 72, 178, 13.03 , 0.21 },
280  { 73, 181, 12.86 , 0.13 },
281  { 74, 184, 12.76 , 0.16 },
282  { 79, 197, 13.35 , 0.10 },
283  { 80, 201, 12.74 , 0.18 },
284  { 81, 205, 13.85 , 0.17 },
285  { 82, 207, 13.295 , 0.071 },
286  { 83, 209, 13.238 , 0.065 },
287  { 90, 232, 12.555 , 0.049 },
288  { 92, 238, 12.592 , 0.035 },
289  { 92, 233, 14.27 , 0.15 },
290  { 92, 235, 13.470 , 0.085 },
291  { 92, 236, 13.90 , 0.40 },
292  { 93, 237, 13.58 , 0.18 },
293  { 94, 239, 13.90 , 0.20 },
294  { 94, 242, 12.86 , 0.19 }
295  };
296 
297  G4double lambda = -1.;
298 
299  size_t nCapRates = sizeof(capRates)/sizeof(capRates[0]);
300  for (size_t j = 0; j < nCapRates; ++j) {
301  if( capRates[j].Z == Z && capRates[j].A == A ) {
302  lambda = capRates[j].cRate / microsecond;
303  break;
304  }
305  // make sure the data is sorted for the next statement to work correctly
306  if (capRates[j].Z > Z) {break;}
307  }
308 
309  if (lambda < 0.) {
310 
311  // == Mu capture lifetime (Goulard and Primakoff PRC10(1974)2034.
312 
313  static const G4double b0a = -0.03;
314  static const G4double b0b = -0.25;
315  static const G4double b0c = 3.24;
316  static const G4double t1 = 875.e-9; // -10-> -9 suggested by user
317  G4double r1 = GetMuonZeff(Z);
318  G4double zeff2 = r1 * r1;
319 
320  // ^-4 -> ^-5 suggested by user
321  G4double xmu = zeff2 * 2.663e-5;
322  G4double a2ze = 0.5 *G4double(A) / G4double(Z);
323  G4double r2 = 1.0 - xmu;
324  lambda = t1 * zeff2 * zeff2 * (r2 * r2) * (1.0 - (1.0 - xmu) * .75704) *
325  (a2ze * b0a + 1.0 - (a2ze - 1.0) * b0b -
326  G4double(2 * (A - Z) + std::abs(a2ze - 1.) ) * b0c / G4double(A * 4) );
327 
328  }
329 
330  return lambda;
331 
332 }
333 
334 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
335 
336 
337 G4double G4MuonMinusBoundDecay::GetMuonZeff(G4int ZZ)
338 {
339 
340  // == Effective charges from
341  // "Total Nuclear Capture Rates for Negative Muons"
342  // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
343  // and if not present from
344  // Ford and Wills Nucl Phys 35(1962)295 or interpolated
345 
346  static const G4int maxZ = 100;
347  static const G4double zeff[] =
348  { 0.,
349  1.00, 1.98, 2.94, 3.89, 4.81, 5.72, 6.61, 7.49, 8.32, 9.14,
350  9.95,10.69,11.48,12.22,12.90,13.64,14.24,14.89,15.53,16.15,
351  16.77,17.38,18.04,18.49,19.06,19.59,20.13,20.66,21.12,21.61,
352  22.02,22.43,22.84,23.24,23.65,24.06,24.47,24.85,25.23,25.61,
353  25.99,26.37,26.69,27.00,27.32,27.63,27.95,28.20,28.42,28.64,
354  28.79,29.03,29.27,29.51,29.75,29.99,30.22,30.36,30.53,30.69,
355  30.85,31.01,31.18,31.34,31.48,31.62,31.76,31.90,32.05,32.19,
356  32.33,32.47,32.61,32.76,32.94,33.11,33.29,33.46,33.64,33.81,
357  34.21,34.18,34.00,34.10,34.21,34.31,34.42,34.52,34.63,34.73,
358  34.84,34.94,35.05,35.16,35.25,35.36,35.46,35.57,35.67,35.78 };
359 
360  G4int Z = std::max(std::min(ZZ, maxZ), 1);
361  return zeff[Z];
362 }
363 
364 G4double G4MuonMinusBoundDecay::GetMuonDecayRate(G4int Z)
365 {
366  G4int A = G4lrint(G4NistManager::Instance()->GetAtomicMassAmu(Z));
367  return GetMuonDecayRate(Z, A, G4MuonMinus::MuonMinus()->GetPDGMass(),
368  G4NucleiProperties::GetNuclearMass(A,Z));
369 }
370 
371 G4double G4MuonMinusBoundDecay::GetMuonDecayRate(G4int Z, G4int A,
372  G4double muMass,
373  G4double nuclMass)
374 {
375  // Decay time correction based on
376  // H. C. Von Baeyer and D. Leiter, Phys. Rev. A 19, 1371 (1979).
377  // replacing N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1.
378  // for Z < 14 inspired by report 2049
379  // Lambda(bound)/Lambda(free) = 1-(0.5+0.06*Mu_Mass/NuclMass)(Z*alpha)^2
380 
381  // PDG 2012 muon lifetime value is 2.1969811(22) 1e-6s
382  // which when inverted gives 0.45517005 1e+6/s
383 
384  // we pass known alraedy in ApplyYourself muon and nucleus masses to
385  // avoid refetching/recalculations
386 
387  struct decRate {
388  G4int Z;
389  G4int A;
390  G4double dRate;
391  G4double dRErr;
392  };
393 
394  // this struct has to be sorted by Z when initialized as we exit the
395  // loop once Z is above the stored value
396 
397  static const decRate decRates [] = {
398  { 0, 0, 0.45517005, 0.00000046 } // free muon
399  };
400 
401  G4double lambda = -1.;
402  if (Z == 0 && A == 0) {lambda = decRates[0].dRate/microsecond;} // free muon
403 
404  // size_t nDecRates = sizeof(decRates)/sizeof(decRates[0]);
405  // for (size_t j = 0; j < nDecRates; ++j) {
406  // if( decRates[j].Z == Z && decRates[j].A == A ) {
407  // lambda = decRates[j].cRate / microsecond;
408  // break;
409  // }
410  // // make sure the data is sorted for the next statement to work correctly
411  // if (decRates[j].Z > Z) {break;}
412  // }
413 
414  // we'll use the above code once we have the data
415  // if we had one value we would just assign it
416  // if (Z == 1 && A == 1) {lambda = decRates[1].dRate/microsecond;}
417 
418  if (lambda < 0.) {
419  const G4double freeMuonDecayRate = 0.45517005 / microsecond;
420  lambda = 1.0;
421  G4double x = Z*fine_structure_const;
422  if (Z<14) {
423  // using the Phys. Rev. formula
424  lambda -= x * x * (0.5 + 0.06 * muMass/nuclMass);
425  } else {
426  // based on a fit to the data ref. in Phys.Rev. C35 (1987) 2212
427  lambda -= x * x * (0.868699 - x * 0.708985);
428  }
429  lambda *= freeMuonDecayRate;
430  }
431  return lambda;
432 }
433 
434 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
435 
436 void G4MuonMinusBoundDecay::ModelDescription(std::ostream& outFile) const
437 {
438  outFile << " Sample probabilities of mu- nuclear capture of decay"
439  << " from K-shell orbit.\n"
440  << " Time of projectile is changed taking into account life time"
441  << " of muonic atom.\n"
442  << " If decay is sampled primary state become stopAndKill,"
443  << " else - isAlive.\n"
444  << " Mainly based on:\n"
445  << " H.C. Von Baeyer and D.Leiter, Phys. Rev. A 19, 1371 (1979)\n"
446  << " T.Suzuki, D.F.Measday, J.P.Roalsvig Phys.Rev. C35 (1987) 2212\n"
447  << " with an emprical fit to the Huff factors for Z >= 14\n"
448  << " from the above review\n";
449 }
450 
451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....