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.1)


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