Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/stopping/src/G4MuonMinusBoundDecay.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/stopping/src/G4MuonMinusBoundDecay.cc (Version 11.3.0) and /processes/hadronic/stopping/src/G4MuonMinusBoundDecay.cc (Version 10.4.p2)


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