Geant4: G4MuonMinusBoundDecay Class Reference

G4MuonMinusBoundDecay::G4MuonMinusBoundDecay ( )

Definition at line 62 of file G4MuonMinusBoundDecay.cc.

References G4ParticleDefinition::GetPDGMass(), and G4MuonMinus::MuonMinus().

00063   : G4HadronicInteraction("muMinusBoundDeacy")
00064 {
00065   fMuMass = G4MuonMinus::MuonMinus()->GetPDGMass(); 
00066 }

G4MuonMinusBoundDecay::~G4MuonMinusBoundDecay ( )

Definition at line 70 of file G4MuonMinusBoundDecay.cc.

00071 {}

G4HadFinalState * G4MuonMinusBoundDecay::ApplyYourself	(	const G4HadProjectile &	aTrack,
		G4Nucleus &	targetNucleus
	)			`[virtual]`

Implements G4HadronicInteraction.

Definition at line 76 of file G4MuonMinusBoundDecay.cc.

References G4AntiNeutrinoE::AntiNeutrinoE(), G4HadFinalState::Clear(), G4Electron::Electron(), G4RandomDirection(), G4UniformRand, G4Nucleus::GetA_asInt(), G4HadProjectile::GetBoundEnergy(), GetMuonCaptureRate(), GetMuonDecayRate(), G4Nucleus::GetZ_asInt(), isAlive, G4InuclParticleNames::lambda, G4NeutrinoMu::NeutrinoMu(), G4HadProjectile::SetGlobalTime(), G4HadFinalState::SetStatusChange(), and stopAndKill.

00078 {
00079   result.Clear();
00080   G4int Z = targetNucleus.GetZ_asInt(); 
00081   G4int A = targetNucleus.GetA_asInt(); 
00082 
00083   // Decide on Decay or Capture, and doit.
00084   G4double lambdac  = GetMuonCaptureRate(Z, A);
00085   G4double lambdad  = GetMuonDecayRate(Z);
00086   G4double lambda   = lambdac + lambdad;
00087 
00088   // ===  sample capture  time and change time of projectile
00089 
00090   G4double time = -std::log(G4UniformRand()) / lambda;
00091   G4HadProjectile* p = const_cast<G4HadProjectile*>(&projectile);
00092   p->SetGlobalTime(time);
00093     
00094   //G4cout << "lambda= " << lambda << " lambdac= " << lambdac 
00095   //<< " t= " << time << G4endl;
00096  
00097   // cascade
00098   if( G4UniformRand()*lambda < lambdac) {
00099     result.SetStatusChange(isAlive);
00100 
00101   } else {
00102 
00103     // Simulation on Decay of mu- on a K-shell of the muonic atom
00104     result.SetStatusChange(stopAndKill);
00105     G4double xmax = 1 + electron_mass_c2*electron_mass_c2/(fMuMass*fMuMass);
00106     G4double xmin = 2.0*electron_mass_c2/fMuMass;
00107     G4double KEnergy = projectile.GetBoundEnergy();
00108 
00109     /*
00110       G4cout << "G4MuonMinusBoundDecay::ApplyYourself" 
00111       << " XMAX= " << xmax << " Ebound= " << KEnergy<< G4endl;
00112     */
00113     G4double pmu = std::sqrt(KEnergy*(KEnergy + 2.0*fMuMass));
00114     G4double emu = KEnergy + fMuMass;
00115     G4ThreeVector dir = G4RandomDirection();
00116     G4LorentzVector MU(pmu*dir, emu);
00117     G4ThreeVector bst = MU.boostVector();
00118 
00119     G4double Eelect, Pelect, x, ecm;
00120     G4LorentzVector EL, NN;
00121     // Calculate electron energy
00122     do {
00123       do {
00124         x = xmin + (xmax-xmin)*G4UniformRand();
00125       } while (G4UniformRand() > (3.0 - 2.0*x)*x*x );
00126       Eelect = x*fMuMass*0.5;
00127       Pelect = 0.0;
00128       if(Eelect > electron_mass_c2) { 
00129         Pelect = std::sqrt(Eelect*Eelect - electron_mass_c2*electron_mass_c2);
00130       } else {
00131         Pelect = 0.0;
00132         Eelect = electron_mass_c2;
00133       }
00134       dir = G4RandomDirection();
00135       EL = G4LorentzVector(Pelect*dir,Eelect);
00136       EL.boost(bst);
00137       Eelect = EL.e() - electron_mass_c2 - 2.0*KEnergy;
00138       //
00139       // Calculate rest frame parameters of 2 neutrinos
00140       //
00141       NN = MU - EL;
00142       ecm = NN.mag2();
00143     } while (Eelect < 0.0 || ecm < 0.0);
00144 
00145     //
00146     // Create electron
00147     //
00148     G4DynamicParticle* dp = new G4DynamicParticle(G4Electron::Electron(),
00149                                                   EL.vect().unit(),
00150                                                   Eelect);
00151 
00152     AddNewParticle(dp, time);
00153     //
00154     // Create Neutrinos
00155     //
00156     ecm = 0.5*std::sqrt(ecm);
00157     bst = NN.boostVector();
00158     G4ThreeVector p1 = ecm * G4RandomDirection();
00159     G4LorentzVector N1 = G4LorentzVector(p1,ecm);
00160     N1.boost(bst);
00161     dp = new G4DynamicParticle(G4AntiNeutrinoE::AntiNeutrinoE(), N1);
00162     AddNewParticle(dp, time);
00163     NN -= N1;
00164     dp = new G4DynamicParticle(G4NeutrinoMu::NeutrinoMu(), NN);
00165     AddNewParticle(dp, time);
00166   }
00167   return &result;
00168 }

G4double G4MuonMinusBoundDecay::GetMuonCaptureRate	(	G4int	Z,
		G4int	A
	)			`[static]`

Definition at line 172 of file G4MuonMinusBoundDecay.cc.

References GetMuonZeff(), and G4InuclParticleNames::lambda.

Referenced by ApplyYourself(), and G4StopElementSelector::GetMuonCaptureRate().

00173 {
00174 
00175   // Initialize data
00176 
00177   // Mu- capture data from 
00178   // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
00179   // weighted average of the two most precise measurements
00180 
00181   // Data for Hydrogen from Phys. Rev. Lett. 99(2007)032002
00182   // Data for Helium from D.F. Measday Phys. Rep. 354(2001)243
00183 
00184   struct capRate {
00185     G4int        Z;
00186     G4int        A;
00187     G4double cRate;
00188     G4double cRErr;
00189   };
00190 
00191   // this struct has to be sorted by Z when initialized as we exit the
00192   // loop once Z is above the stored value; cRErr are not used now but
00193   // are included for completeness and future use
00194 
00195   const capRate capRates [] = {
00196     {  1,   1,  0.000725, 0.000017 },
00197     {  2,   3,  0.002149, 0.00017 }, 
00198     {  2,   4,  0.000356, 0.000026 },
00199     {  3,   6,  0.004647, 0.00012 },
00200     {  3,   7,  0.002229, 0.00012 },
00201     {  4,   9,  0.006107, 0.00019 },
00202     {  5,  10,  0.02757 , 0.00063 },
00203     {  5,  11,  0.02188 , 0.00064 },
00204     {  6,  12,  0.03807 , 0.00031 },
00205     {  6,  13,  0.03474 , 0.00034 },
00206     {  7,  14,  0.06885 , 0.00057 },
00207     {  8,  16,  0.10242 , 0.00059 },
00208     {  8,  18,  0.0880  , 0.0015  },
00209     {  9,  19,  0.22905 , 0.00099 },
00210     { 10,  20,  0.2288  , 0.0045 },
00211     { 11,  23,  0.3773  , 0.0014 },
00212     { 12,  24,  0.4823  , 0.0013 },
00213     { 13,  27,  0.6985  , 0.0012 },
00214     { 14,  28,  0.8656  , 0.0015 },
00215     { 15,  31,  1.1681  , 0.0026 },
00216     { 16,  32,  1.3510  , 0.0029 },
00217     { 17,  35,  1.800   , 0.050 },
00218     { 17,  37,  1.250   , 0.050 },
00219     { 18,  40,  1.2727  , 0.0650 },
00220     { 19,  39,  1.8492  , 0.0050 },
00221     { 20,  40,  2.5359  , 0.0070 },
00222     { 21,  45,  2.711   , 0.025 },
00223     { 22,  48,  2.5908  , 0.0115 },
00224     { 23,  51,  3.073   , 0.022 },
00225     { 24,  50,  3.825   , 0.050 },
00226     { 24,  52,  3.465   , 0.026 },
00227     { 24,  53,  3.297   , 0.045 },
00228     { 24,  54,  3.057   , 0.042 },
00229     { 25,  55,  3.900   , 0.030 },
00230     { 26,  56,  4.408   , 0.022 },
00231     { 27,  59,  4.945   , 0.025 },
00232     { 28,  58,  6.11    , 0.10 },
00233     { 28,  60,  5.56    , 0.10 },
00234     { 28,  62,  4.72    , 0.10 },
00235     { 29,  63,  5.691   , 0.030 },
00236     { 30,  66,  5.806   , 0.031 },
00237     { 31,  69,  5.700   , 0.060 },
00238     { 32,  72,  5.561   , 0.031 },
00239     { 33,  75,  6.094   , 0.037 },
00240     { 34,  80,  5.687   , 0.030 },
00241     { 35,  79,  7.223   , 0.28 },
00242     { 35,  81,  7.547   , 0.48 },
00243     { 37,  85,  6.89    , 0.14 },
00244     { 38,  88,  6.93    , 0.12 },
00245     { 39,  89,  7.89    , 0.11 },
00246     { 40,  91,  8.620   , 0.053 },
00247     { 41,  93, 10.38    , 0.11 },
00248     { 42,  96,  9.298   , 0.063 },
00249     { 45, 103, 10.010   , 0.045 },
00250     { 46, 106, 10.000   , 0.070 },
00251     { 47, 107, 10.869   , 0.095 },
00252     { 48, 112, 10.624   , 0.094 },
00253     { 49, 115, 11.38    , 0.11 },
00254     { 50, 119, 10.60    , 0.11 },
00255     { 51, 121, 10.40    , 0.12 },
00256     { 52, 128,  9.174   , 0.074 },
00257     { 53, 127, 11.276   , 0.098 },
00258     { 55, 133, 10.98    , 0.25 },
00259     { 56, 138, 10.112   , 0.085 },
00260     { 57, 139, 10.71    , 0.10 },
00261     { 58, 140, 11.501   , 0.087 },
00262     { 59, 141, 13.45    , 0.13 },
00263     { 60, 144, 12.35    , 0.13 },
00264     { 62, 150, 12.22    , 0.17 },
00265     { 64, 157, 12.00    , 0.13 },
00266     { 65, 159, 12.73    , 0.13 },
00267     { 66, 163, 12.29    , 0.18 },
00268     { 67, 165, 12.95    , 0.13 },
00269     { 68, 167, 13.04    , 0.27 },
00270     { 72, 178, 13.03    , 0.21 },
00271     { 73, 181, 12.86    , 0.13 },
00272     { 74, 184, 12.76    , 0.16 },
00273     { 79, 197, 13.35    , 0.10 },
00274     { 80, 201, 12.74    , 0.18 },
00275     { 81, 205, 13.85    , 0.17 },
00276     { 82, 207, 13.295   , 0.071 },
00277     { 83, 209, 13.238   , 0.065 },
00278     { 90, 232, 12.555   , 0.049 },
00279     { 92, 238, 12.592   , 0.035 },
00280     { 92, 233, 14.27    , 0.15 },
00281     { 92, 235, 13.470   , 0.085 },
00282     { 92, 236, 13.90    , 0.40 },
00283     { 93, 237, 13.58    , 0.18 },
00284     { 94, 239, 13.90    , 0.20 },
00285     { 94, 242, 12.86    , 0.19 }
00286   };
00287 
00288   G4double lambda = -1.;
00289 
00290   size_t nCapRates = sizeof(capRates)/sizeof(capRates[0]);
00291   for (size_t j = 0; j < nCapRates; ++j) {
00292     if( capRates[j].Z == Z && capRates[j].A == A ) {
00293       lambda = capRates[j].cRate / microsecond;
00294       break;
00295     }
00296     // make sure the data is sorted for the next statement to work correctly
00297     if (capRates[j].Z > Z) {break;}
00298   }
00299 
00300   if (lambda < 0.) {
00301 
00302     // ==  Mu capture lifetime (Goulard and Primakoff PRC10(1974)2034.
00303 
00304     const G4double b0a = -0.03;
00305     const G4double b0b = -0.25;
00306     const G4double b0c =  3.24;
00307     const G4double t1 = 875.e-9; // -10-> -9  suggested by user
00308     G4double r1 = GetMuonZeff(Z);
00309     G4double zeff2 = r1 * r1;
00310 
00311     // ^-4 -> ^-5 suggested by user
00312     G4double xmu = zeff2 * 2.663e-5;
00313     G4double a2ze = 0.5 *G4double(A) / G4double(Z);
00314     G4double r2 = 1.0 - xmu;
00315     lambda = t1 * zeff2 * zeff2 * (r2 * r2) * (1.0 - (1.0 - xmu) * .75704) *
00316       (a2ze * b0a + 1.0 - (a2ze - 1.0) * b0b -
00317        G4double(2 * (A - Z)  + std::fabs(a2ze - 1.) ) * b0c / G4double(A * 4) );
00318 
00319   }
00320 
00321   return lambda;
00322 
00323 }

G4double G4MuonMinusBoundDecay::GetMuonDecayRate ( G4int Z ) [static]

Definition at line 360 of file G4MuonMinusBoundDecay.cc.

References GetMuonZeff(), and G4InuclParticleNames::lambda.

Referenced by ApplyYourself(), and G4StopElementSelector::GetMuonDecayRate().

00361 {
00362   // Decay time on K-shell 
00363   // N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1.
00364 
00365   // this is the "small Z" approximation formula (2.9)
00366   // Lambda(bound)/Lambda(free) = 1-beta(Z*alpha)**2 with beta~=2.5
00367   // we assume that Z is Zeff
00368 
00369   // PDG 2012 muon lifetime value is 2.1969811(22) 10e-6s
00370   // which when inverted gives       0.45517005    10e+6/s
00371 
00372   struct decRate {
00373     G4int         Z;
00374     G4double  dRate;
00375     G4double  dRErr;
00376   };
00377 
00378   // this struct has to be sorted by Z when initialized as we exit the
00379   // loop once Z is above the stored value
00380 
00381   const decRate decRates [] = {
00382     {  1,  0.4558514, 0.0000151 }
00383   };
00384 
00385   G4double lambda = -1.;
00386 
00387   // size_t nDecRates = sizeof(decRates)/sizeof(decRates[0]);
00388   // for (size_t j = 0; j < nDecRates; ++j) {
00389   //   if( decRates[j].Z == Z ) {
00390   //     lambda = decRates[j].dRate / microsecond;
00391   //     break;
00392   //   }
00393   //   // make sure the data is sorted for the next statement to work
00394   //   if (decRates[j].Z > Z) {break;}
00395   // }
00396 
00397   // we'll use the above code once we have more data
00398   // since we only have one value we just assign it
00399   if (Z == 1) {lambda =  decRates[0].dRate/microsecond;}
00400 
00401   if (lambda < 0.) {
00402     const G4double freeMuonDecayRate =  0.45517005 / microsecond;
00403     lambda = 1.0;
00404     G4double x = GetMuonZeff(Z)*fine_structure_const;
00405     lambda -= 2.5 * x * x;
00406     lambda *= freeMuonDecayRate;
00407   }
00408 
00409   return lambda;
00410 
00411 }

G4double G4MuonMinusBoundDecay::GetMuonZeff ( G4int Z ) [static]

Definition at line 328 of file G4MuonMinusBoundDecay.cc.

Referenced by GetMuonCaptureRate(), and GetMuonDecayRate().

00329 {
00330 
00331   // ==  Effective charges from 
00332   // "Total Nuclear Capture Rates for Negative Muons"
00333   // T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212
00334   // and if not present from
00335   // Ford and Wills Nucl Phys 35(1962)295 or interpolated
00336 
00337 
00338   const size_t maxZ = 100;
00339   const G4double zeff[maxZ+1] = 
00340     {  0.,
00341        1.00, 1.98, 2.94, 3.89, 4.81, 5.72, 6.61, 7.49, 8.32, 9.14,
00342        9.95,10.69,11.48,12.22,12.90,13.64,14.24,14.89,15.53,16.15,
00343        16.77,17.38,18.04,18.49,19.06,19.59,20.13,20.66,21.12,21.61,
00344        22.02,22.43,22.84,23.24,23.65,24.06,24.47,24.85,25.23,25.61,
00345        25.99,26.37,26.69,27.00,27.32,27.63,27.95,28.20,28.42,28.64,
00346        28.79,29.03,29.27,29.51,29.75,29.99,30.22,30.36,30.53,30.69,
00347        30.85,31.01,31.18,31.34,31.48,31.62,31.76,31.90,32.05,32.19,
00348        32.33,32.47,32.61,32.76,32.94,33.11,33.29,33.46,33.64,33.81,
00349        34.21,34.18,34.00,34.10,34.21,34.31,34.42,34.52,34.63,34.73,
00350        34.84,34.94,35.05,35.16,35.25,35.36,35.46,35.57,35.67,35.78 };
00351 
00352   if (Z<0) {Z=0;}
00353   if (Z>G4int(maxZ)) {Z=maxZ;}
00354 
00355   return zeff[Z];
00356 
00357 }

void G4MuonMinusBoundDecay::ModelDescription ( std::ostream & outFile ) const [virtual]

Reimplemented from G4HadronicInteraction.

Definition at line 415 of file G4MuonMinusBoundDecay.cc.

00416 {
00417   outFile << "Sample probabilities of mu- nuclear capture of decay"
00418           << " from K-shell orbit.\n"
00419           << "   Time of projectile is changed taking into account life time"
00420           << " of muonic atom.\n"
00421           << "   If decay is sampled primary state become stopAndKill,"
00422           << " else - isAlive.\n"
00423           << "  Based of reviews:\n"
00424           << "  N.C.Mukhopadhyay Phy. Rep. 30 (1977) 1.\n"
00425           << "  T. Suzuki, D. F. Measday, J.P. Roalsvig Phys.Rev. C35 (1987) 2212\n"; 
00426 
00427 }


Public Member Functions
	G4MuonMinusBoundDecay ()
	~G4MuonMinusBoundDecay ()
G4HadFinalState *	ApplyYourself (const G4HadProjectile &aTrack, G4Nucleus &targetNucleus)
void	ModelDescription (std::ostream &outFile) const
Static Public Member Functions
static G4double	GetMuonCaptureRate (G4int Z, G4int A)
static G4double	GetMuonDecayRate (G4int Z)
static G4double	GetMuonZeff (G4int Z)

G4MuonMinusBoundDecay Class Reference

Public Member Functions

Static Public Member Functions

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation