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
00046
00047
00048 #include "G4KleinNishinaModel.hh"
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4Electron.hh"
00052 #include "G4Gamma.hh"
00053 #include "Randomize.hh"
00054 #include "G4RandomDirection.hh"
00055 #include "G4DataVector.hh"
00056 #include "G4ParticleChangeForGamma.hh"
00057 #include "G4VAtomDeexcitation.hh"
00058 #include "G4AtomicShells.hh"
00059 #include "G4LossTableManager.hh"
00060
00061
00062
00063 using namespace std;
00064
00065 G4KleinNishinaModel::G4KleinNishinaModel(const G4String& nam)
00066 : G4VEmModel(nam)
00067 {
00068 theGamma = G4Gamma::Gamma();
00069 theElectron = G4Electron::Electron();
00070 lowestGammaEnergy = 1.0*eV;
00071 limitFactor = 4;
00072 fProbabilities.resize(9,0.0);
00073 SetDeexcitationFlag(true);
00074 fParticleChange = 0;
00075 fAtomDeexcitation = 0;
00076 }
00077
00078
00079
00080 G4KleinNishinaModel::~G4KleinNishinaModel()
00081 {}
00082
00083
00084
00085 void G4KleinNishinaModel::Initialise(const G4ParticleDefinition* p,
00086 const G4DataVector& cuts)
00087 {
00088 fAtomDeexcitation = G4LossTableManager::Instance()->AtomDeexcitation();
00089 InitialiseElementSelectors(p, cuts);
00090 if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
00091 }
00092
00093
00094
00095 G4double
00096 G4KleinNishinaModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00097 G4double GammaEnergy,
00098 G4double Z, G4double,
00099 G4double, G4double)
00100 {
00101 G4double xSection = 0.0 ;
00102 if ( Z < 0.9999 || GammaEnergy < 0.1*keV) { return xSection; }
00103
00104 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
00105
00106 static const G4double
00107 d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527 *barn, d4=-1.9798e+1*barn,
00108 e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
00109 f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
00110
00111 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
00112 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
00113
00114 G4double T0 = 15.0*keV;
00115 if (Z < 1.5) { T0 = 40.0*keV; }
00116
00117 G4double X = max(GammaEnergy, T0) / electron_mass_c2;
00118 xSection = p1Z*std::log(1.+2.*X)/X
00119 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
00120
00121
00122 if (GammaEnergy < T0) {
00123 G4double dT0 = keV;
00124 X = (T0+dT0) / electron_mass_c2 ;
00125 G4double sigma = p1Z*log(1.+2*X)/X
00126 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
00127 G4double c1 = -T0*(sigma-xSection)/(xSection*dT0);
00128 G4double c2 = 0.150;
00129 if (Z > 1.5) { c2 = 0.375-0.0556*log(Z); }
00130 G4double y = log(GammaEnergy/T0);
00131 xSection *= exp(-y*(c1+c2*y));
00132 }
00133
00134 if(xSection < 0.0) { xSection = 0.0; }
00135
00136
00137 return xSection;
00138 }
00139
00140
00141
00142 void G4KleinNishinaModel::SampleSecondaries(
00143 std::vector<G4DynamicParticle*>* fvect,
00144 const G4MaterialCutsCouple* couple,
00145 const G4DynamicParticle* aDynamicGamma,
00146 G4double,
00147 G4double)
00148 {
00149
00150 G4double energy = aDynamicGamma->GetKineticEnergy();
00151 G4ThreeVector direction = aDynamicGamma->GetMomentumDirection();
00152
00153
00154 const G4Element* elm = SelectRandomAtom(couple, theGamma, energy);
00155
00156
00157 G4int nShells = elm->GetNbOfAtomicShells();
00158 if(nShells > (G4int)fProbabilities.size()) { fProbabilities.resize(nShells); }
00159 G4double totprob = 0.0;
00160 G4int i;
00161 for(i=0; i<nShells; ++i) {
00162
00163 totprob += elm->GetNbOfShellElectrons(i);
00164
00165 fProbabilities[i] = totprob;
00166 }
00167
00168
00169
00170
00171
00172
00173 G4double bindingEnergy, ePotEnergy, eKinEnergy;
00174 G4double gamEnergy0, gamEnergy1;
00175
00176
00177
00178 do {
00179
00180 G4double xprob = totprob*G4UniformRand();
00181
00182
00183 for(i=0; i<nShells; ++i) { if(xprob <= fProbabilities[i]) { break; } }
00184
00185 bindingEnergy = elm->GetAtomicShell(i);
00186
00187
00188 lv1.set(0.0,0.0,energy,energy);
00189
00190
00191
00192
00193
00194
00195
00196 G4double x = -log(G4UniformRand());
00197 eKinEnergy = bindingEnergy*x;
00198 ePotEnergy = bindingEnergy*(1.0 + x);
00199
00200
00201 G4double eTotMomentum = sqrt(eKinEnergy*(eKinEnergy + 2*electron_mass_c2));
00202 G4double phi = G4UniformRand()*twopi;
00203 G4double costet = 2*G4UniformRand() - 1;
00204 G4double sintet = sqrt((1 - costet)*(1 + costet));
00205 lv2.set(eTotMomentum*sintet*cos(phi),eTotMomentum*sintet*sin(phi),
00206 eTotMomentum*costet,eKinEnergy + electron_mass_c2);
00207 bst = lv2.boostVector();
00208 lv1.boost(-bst);
00209
00210 gamEnergy0 = lv1.e();
00211
00212
00213
00214
00215
00216 G4double E0_m = gamEnergy0/electron_mass_c2;
00217
00218
00219
00220
00221
00222 G4double epsilon, epsilonsq, onecost, sint2, greject ;
00223
00224 G4double eps0 = 1./(1 + 2*E0_m);
00225 G4double epsilon0sq = eps0*eps0;
00226 G4double alpha1 = - log(eps0);
00227 G4double alpha2 = 0.5*(1 - epsilon0sq);
00228
00229 do {
00230 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) {
00231 epsilon = exp(-alpha1*G4UniformRand());
00232 epsilonsq = epsilon*epsilon;
00233
00234 } else {
00235 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
00236 epsilon = sqrt(epsilonsq);
00237 }
00238
00239 onecost = (1.- epsilon)/(epsilon*E0_m);
00240 sint2 = onecost*(2.-onecost);
00241 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
00242
00243 } while (greject < G4UniformRand());
00244 gamEnergy1 = epsilon*gamEnergy0;
00245
00246
00247 lv2.set(0.0,0.0,0.0,electron_mass_c2);
00248 lv2 += lv1;
00249
00250
00251
00252
00253 if(sint2 < 0.0) { sint2 = 0.0; }
00254 costet = 1. - onecost;
00255 sintet = sqrt(sint2);
00256 phi = twopi * G4UniformRand();
00257
00258
00259
00260
00261 G4ThreeVector gamDir = lv1.vect().unit();
00262 G4ThreeVector v = G4ThreeVector(sintet*cos(phi),sintet*sin(phi),costet);
00263 v.rotateUz(gamDir);
00264 lv1.set(gamEnergy1*v.x(),gamEnergy1*v.y(),gamEnergy1*v.z(),gamEnergy1);
00265 lv2 -= lv1;
00266
00267 lv2.boost(bst);
00268 eKinEnergy = lv2.e() - electron_mass_c2 - ePotEnergy;
00269
00270
00271 } while ( eKinEnergy < 0.0 );
00272
00273
00274
00275
00276
00277 lv1.boost(bst);
00278 gamEnergy1 = lv1.e();
00279 if(gamEnergy1 > lowestGammaEnergy) {
00280 G4ThreeVector gamDirection1 = lv1.vect().unit();
00281 gamDirection1.rotateUz(direction);
00282 fParticleChange->ProposeMomentumDirection(gamDirection1);
00283 } else {
00284 fParticleChange->ProposeTrackStatus(fStopAndKill);
00285 gamEnergy1 = 0.0;
00286 }
00287 fParticleChange->SetProposedKineticEnergy(gamEnergy1);
00288
00289
00290
00291
00292
00293 if(eKinEnergy > lowestGammaEnergy) {
00294 G4ThreeVector eDirection = lv2.vect().unit();
00295 eDirection.rotateUz(direction);
00296 G4DynamicParticle* dp =
00297 new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
00298 fvect->push_back(dp);
00299 } else { eKinEnergy = 0.0; }
00300
00301 G4double edep = energy - gamEnergy1 - eKinEnergy;
00302
00303
00304
00305 if(fAtomDeexcitation) {
00306 G4int index = couple->GetIndex();
00307 if(fAtomDeexcitation->CheckDeexcitationActiveRegion(index)) {
00308 G4int Z = G4lrint(elm->GetZ());
00309 G4AtomicShellEnumerator as = G4AtomicShellEnumerator(i);
00310 const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
00311 size_t nbefore = fvect->size();
00312 fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
00313 size_t nafter = fvect->size();
00314 if(nafter > nbefore) {
00315 for (size_t j=nbefore; j<nafter; ++j) {
00316 edep -= ((*fvect)[j])->GetKineticEnergy();
00317 }
00318 }
00319 }
00320 }
00321
00322 if(edep < 0.0) { edep = 0.0; }
00323 fParticleChange->ProposeLocalEnergyDeposit(edep);
00324 }
00325
00326
00327