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