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