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
00049 #include <CLHEP/Random/RandGamma.h>
00050 #include "globals.hh"
00051 #include "G4PhysicalConstants.hh"
00052 #include "G4SystemOfUnits.hh"
00053 #include "G4RandomDirection.hh"
00054 #include "Randomize.hh"
00055
00056 #include "G4HeatedKleinNishinaCompton.hh"
00057 #include "G4Electron.hh"
00058 #include "G4Gamma.hh"
00059 #include "Randomize.hh"
00060 #include "G4DataVector.hh"
00061 #include "G4ParticleChangeForGamma.hh"
00062
00063
00064
00065 using namespace std;
00066
00067 G4HeatedKleinNishinaCompton::G4HeatedKleinNishinaCompton(const G4ParticleDefinition*,
00068 const G4String& nam)
00069 : G4VEmModel(nam)
00070 {
00071 theGamma = G4Gamma::Gamma();
00072 theElectron = G4Electron::Electron();
00073 lowestGammaEnergy = 1.0*eV;
00074 fTemperature = 1.0*keV;
00075 fParticleChange = 0;
00076 }
00077
00078
00079
00080 G4HeatedKleinNishinaCompton::~G4HeatedKleinNishinaCompton()
00081 {}
00082
00083
00084
00085 void G4HeatedKleinNishinaCompton::Initialise(const G4ParticleDefinition*,
00086 const G4DataVector&)
00087 {
00088 if(!fParticleChange) fParticleChange = GetParticleChangeForGamma();
00089 }
00090
00092
00093
00094
00095 G4double G4HeatedKleinNishinaCompton::ComputeCrossSectionPerAtom(
00096 const G4ParticleDefinition*,
00097 G4double GammaEnergy,
00098 G4double Z, G4double,
00099 G4double, G4double)
00100 {
00101 G4double xSection = 0.0 ;
00102 if ( Z < 0.9999 ) return xSection;
00103 if ( GammaEnergy < 0.01*keV ) return xSection;
00104
00105
00106 static const G4double a = 20.0 , b = 230.0 , c = 440.0;
00107
00108 static const G4double
00109 d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527 *barn, d4=-1.9798e+1*barn,
00110 e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
00111 f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
00112
00113 G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
00114 p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
00115
00116 G4double T0 = 15.0*keV;
00117 if (Z < 1.5) T0 = 40.0*keV;
00118
00119 G4double X = max(GammaEnergy, T0) / electron_mass_c2;
00120 xSection = p1Z*std::log(1.+2.*X)/X
00121 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
00122
00123
00124 if (GammaEnergy < T0) {
00125 G4double dT0 = 1.*keV;
00126 X = (T0+dT0) / electron_mass_c2 ;
00127 G4double sigma = p1Z*log(1.+2*X)/X
00128 + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
00129 G4double c1 = -T0*(sigma-xSection)/(xSection*dT0);
00130 G4double c2 = 0.150;
00131 if (Z > 1.5) c2 = 0.375-0.0556*log(Z);
00132 G4double y = log(GammaEnergy/T0);
00133 xSection *= exp(-y*(c1+c2*y));
00134 }
00135
00136 return xSection;
00137 }
00138
00140
00141
00142
00143 void G4HeatedKleinNishinaCompton::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00144 const G4MaterialCutsCouple*,
00145 const G4DynamicParticle* aDynamicGamma,
00146 G4double,
00147 G4double)
00148 {
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158 G4double eMomentumC2 = CLHEP::RandGamma::shoot(1.5,1.);
00159 eMomentumC2 *= 2*electron_mass_c2*fTemperature;
00160 G4ThreeVector eMomDir = G4RandomDirection();
00161 eMomDir *= std::sqrt(eMomentumC2);
00162 G4double eEnergy = std::sqrt(eMomentumC2+electron_mass_c2*electron_mass_c2);
00163 G4LorentzVector electron4v = G4LorentzVector(eMomDir,eEnergy);
00164 G4ThreeVector bst = electron4v.boostVector();
00165
00166 G4LorentzVector gamma4v = aDynamicGamma->Get4Momentum();
00167 gamma4v.boost(-bst);
00168
00169 G4ThreeVector gammaMomV = gamma4v.vect();
00170 G4double gamEnergy0 = gammaMomV.mag();
00171
00172
00173
00174 G4double E0_m = gamEnergy0 / electron_mass_c2 ;
00175
00176
00177
00178 G4ThreeVector gamDirection0 = gammaMomV/gamEnergy0;
00179
00180
00181
00182
00183 G4double epsilon, epsilonsq, onecost, sint2, greject ;
00184
00185 G4double eps0 = 1./(1. + 2.*E0_m);
00186 G4double epsilon0sq = eps0*eps0;
00187 G4double alpha1 = - log(eps0);
00188 G4double alpha2 = 0.5*(1.- epsilon0sq);
00189
00190 do
00191 {
00192 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
00193 {
00194 epsilon = exp(-alpha1*G4UniformRand());
00195 epsilonsq = epsilon*epsilon;
00196
00197 }
00198 else
00199 {
00200 epsilonsq = epsilon0sq + (1.- epsilon0sq)*G4UniformRand();
00201 epsilon = sqrt(epsilonsq);
00202 };
00203
00204 onecost = (1.- epsilon)/(epsilon*E0_m);
00205 sint2 = onecost*(2.-onecost);
00206 greject = 1. - epsilon*sint2/(1.+ epsilonsq);
00207
00208 } while (greject < G4UniformRand());
00209
00210
00211
00212
00213
00214 G4double cosTeta = 1. - onecost;
00215 G4double sinTeta = sqrt (sint2);
00216 G4double Phi = twopi * G4UniformRand();
00217 G4double dirx = sinTeta*cos(Phi), diry = sinTeta*sin(Phi), dirz = cosTeta;
00218
00219
00220
00221
00222
00223 G4ThreeVector gamDirection1 ( dirx,diry,dirz );
00224 gamDirection1.rotateUz(gamDirection0);
00225 G4double gamEnergy1 = epsilon*gamEnergy0;
00226 gamDirection1 *= gamEnergy1;
00227
00228 G4LorentzVector gamma4vfinal = G4LorentzVector(gamDirection1,gamEnergy1);
00229
00230
00231
00232
00233
00234 G4double eKinEnergy = gamEnergy0 - gamEnergy1;
00235 G4ThreeVector eDirection = gamEnergy0*gamDirection0 - gamEnergy1*gamDirection1;
00236 eDirection = eDirection.unit();
00237 G4double eFinalMom = std::sqrt(eKinEnergy*(eKinEnergy+2*electron_mass_c2));
00238 eDirection *= eFinalMom;
00239 G4LorentzVector e4vfinal = G4LorentzVector(eDirection,gamEnergy1+electron_mass_c2);
00240
00241 gamma4vfinal.boost(bst);
00242 e4vfinal.boost(bst);
00243
00244 gamDirection1 = gamma4vfinal.vect();
00245 gamEnergy1 = gamDirection1.mag();
00246 gamDirection1 /= gamEnergy1;
00247
00248
00249
00250
00251 fParticleChange->SetProposedKineticEnergy(gamEnergy1);
00252
00253 if( gamEnergy1 > lowestGammaEnergy )
00254 {
00255 gamDirection1 /= gamEnergy1;
00256 fParticleChange->ProposeMomentumDirection(gamDirection1);
00257 }
00258 else
00259 {
00260 fParticleChange->ProposeTrackStatus(fStopAndKill);
00261 gamEnergy1 += fParticleChange->GetLocalEnergyDeposit();
00262 fParticleChange->ProposeLocalEnergyDeposit(gamEnergy1);
00263 }
00264
00265 eKinEnergy = e4vfinal.t()-electron_mass_c2;
00266
00267 if( eKinEnergy > DBL_MIN )
00268 {
00269
00270 eDirection = e4vfinal.vect();
00271 G4double eFinMomMag = eDirection.mag();
00272 eDirection /= eFinMomMag;
00273 G4DynamicParticle* dp = new G4DynamicParticle(theElectron,eDirection,eKinEnergy);
00274 fvect->push_back(dp);
00275 }
00276 }
00277
00279
00280