00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045 #include "G4EmCaptureCascade.hh"
00046 #include "G4PhysicalConstants.hh"
00047 #include "G4SystemOfUnits.hh"
00048 #include "Randomize.hh"
00049 #include "G4MuonMinus.hh"
00050 #include "G4Electron.hh"
00051 #include "G4Gamma.hh"
00052 #include "G4NucleiProperties.hh"
00053
00054
00055
00056 G4EmCaptureCascade::G4EmCaptureCascade()
00057 : G4HadronicInteraction("emCaptureCascade")
00058 {
00059 theElectron = G4Electron::Electron();
00060 theGamma = G4Gamma::Gamma();
00061 fMuMass = G4MuonMinus::MuonMinus()->GetPDGMass();
00062 fTime = 0.0;
00063
00064
00065
00066
00067 const G4int nlevels = 28;
00068 const G4int listK[nlevels] = {
00069 1, 2, 4, 6, 8, 11, 14, 17, 18, 21, 24,
00070 26, 29, 32, 38, 40, 41, 44, 49, 53, 55,
00071 60, 65, 70, 75, 81, 85, 92};
00072 const G4double listKEnergy[nlevels] = {
00073 0.00275, 0.011, 0.043, 0.098, 0.173, 0.326,
00074 0.524, 0.765, 0.853, 1.146, 1.472,
00075 1.708, 2.081, 2.475, 3.323, 3.627,
00076 3.779, 4.237, 5.016, 5.647, 5.966,
00077 6.793, 7.602, 8.421, 9.249, 10.222,
00078 10.923,11.984};
00079
00080 fKLevelEnergy[0] = 0.0;
00081 fKLevelEnergy[1] = listKEnergy[0];
00082 G4int idx = 1;
00083 for(G4int i=1; i<nlevels; ++i) {
00084 G4int z1 = listK[idx];
00085 G4int z2 = listK[i];
00086 if(z1+1 < z2) {
00087 G4double dz = G4double(z2 - z1);
00088 G4double y1 = listKEnergy[idx]/G4double(z1*z1);
00089 G4double y2 = listKEnergy[i]/G4double(z2*z2);
00090 for(G4int z=z1+1; z<z2; ++z) {
00091 fKLevelEnergy[z] = (y1 + (y2 - y1)*(z - z1)/dz)*z*z;
00092 }
00093 }
00094 fKLevelEnergy[z2] = listKEnergy[i];
00095 idx = i;
00096 }
00097 for( G4int i = 0; i<14; ++i) { fLevelEnergy[i] = 0.0; }
00098 }
00099
00100
00101
00102 G4EmCaptureCascade::~G4EmCaptureCascade()
00103 {}
00104
00105
00106
00107 G4HadFinalState*
00108 G4EmCaptureCascade::ApplyYourself(const G4HadProjectile& projectile,
00109 G4Nucleus& targetNucleus)
00110 {
00111 result.Clear();
00112 result.SetStatusChange(isAlive);
00113 fTime = projectile.GetGlobalTime();
00114
00115 G4int Z = targetNucleus.GetZ_asInt();
00116 G4int A = targetNucleus.GetA_asInt();
00117 G4double massA = G4NucleiProperties::GetNuclearMass(A, Z);
00118 G4double mass = fMuMass * massA / (fMuMass + massA) ;
00119 G4double e = 13.6 * eV * Z * Z * mass/ electron_mass_c2;
00120
00121
00122 fLevelEnergy[0] = fKLevelEnergy[Z];
00123 for( G4int i = 2; i < 15; ++i) {
00124 fLevelEnergy[i-1] = e/G4double(i*i);
00125 }
00126
00127 G4int nElec = G4int(Z);
00128 G4int nAuger = 1;
00129 G4int nLevel = 13;
00130 G4double pGamma = Z*Z*Z*Z;
00131
00132
00133 G4double edep = fLevelEnergy[13];
00134 AddNewParticle(theElectron,edep);
00135 G4double deltaE;
00136
00137
00138
00139
00140 do {
00141
00142
00143 if((nAuger < nElec) && ((pGamma + 10000.0) * G4UniformRand() < 10000.0) ) {
00144 ++nAuger;
00145 deltaE = fLevelEnergy[nLevel-1] - fLevelEnergy[nLevel];
00146 --nLevel;
00147 AddNewParticle(theElectron, deltaE);
00148
00149 } else {
00150
00151
00152
00153
00154 G4double var = (10.0 + G4double(nLevel - 1) ) * G4UniformRand();
00155 G4int iLevel = nLevel - 1 ;
00156 if(var > 10.0) iLevel -= G4int(var-10.0) + 1;
00157 if( iLevel < 0 ) iLevel = 0;
00158 deltaE = fLevelEnergy[iLevel] - fLevelEnergy[nLevel];
00159 nLevel = iLevel;
00160 AddNewParticle(theGamma, deltaE);
00161 }
00162 edep += deltaE;
00163
00164 } while( nLevel > 0 );
00165
00166 result.SetLocalEnergyDeposit(edep);
00167 return &result;
00168 }
00169
00170
00171
00172 void G4EmCaptureCascade::ModelDescription(std::ostream& outFile) const
00173 {
00174 outFile << "Simulation of electromagnetic cascade from capture level"
00175 << " to K-shell of the mesonic atom\n."
00176 << "Probabilities of gamma and Auger transitions from\n"
00177 << " N.C.Mukhopadhyay Phys. Rep. 30 (1977) 1.\n";
00178 }
00179
00180