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
00050
00051
00052
00053
00054
00055 #include "G4BetheHeitlerModel.hh"
00056 #include "G4PhysicalConstants.hh"
00057 #include "G4SystemOfUnits.hh"
00058 #include "G4Electron.hh"
00059 #include "G4Positron.hh"
00060 #include "G4Gamma.hh"
00061 #include "Randomize.hh"
00062 #include "G4ParticleChangeForGamma.hh"
00063
00064
00065
00066 using namespace std;
00067
00068 G4BetheHeitlerModel::G4BetheHeitlerModel(const G4ParticleDefinition*,
00069 const G4String& nam)
00070 : G4VEmModel(nam)
00071 {
00072 fParticleChange = 0;
00073 theGamma = G4Gamma::Gamma();
00074 thePositron = G4Positron::Positron();
00075 theElectron = G4Electron::Electron();
00076 }
00077
00078
00079
00080 G4BetheHeitlerModel::~G4BetheHeitlerModel()
00081 {}
00082
00083
00084
00085 void G4BetheHeitlerModel::Initialise(const G4ParticleDefinition* p,
00086 const G4DataVector& cuts)
00087 {
00088 if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); }
00089 InitialiseElementSelectors(p, cuts);
00090 }
00091
00092
00093
00094 G4double
00095 G4BetheHeitlerModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00096 G4double GammaEnergy, G4double Z,
00097 G4double, G4double, G4double)
00098
00099
00100
00101
00102
00103
00104 {
00105 static const G4double GammaEnergyLimit = 1.5*MeV;
00106 G4double xSection = 0.0 ;
00107 if ( Z < 0.9 || GammaEnergy <= 2.0*electron_mass_c2 ) { return xSection; }
00108
00109 static const G4double
00110 a0= 8.7842e+2*microbarn, a1=-1.9625e+3*microbarn, a2= 1.2949e+3*microbarn,
00111 a3=-2.0028e+2*microbarn, a4= 1.2575e+1*microbarn, a5=-2.8333e-1*microbarn;
00112
00113 static const G4double
00114 b0=-1.0342e+1*microbarn, b1= 1.7692e+1*microbarn, b2=-8.2381 *microbarn,
00115 b3= 1.3063 *microbarn, b4=-9.0815e-2*microbarn, b5= 2.3586e-3*microbarn;
00116
00117 static const G4double
00118 c0=-4.5263e+2*microbarn, c1= 1.1161e+3*microbarn, c2=-8.6749e+2*microbarn,
00119 c3= 2.1773e+2*microbarn, c4=-2.0467e+1*microbarn, c5= 6.5372e-1*microbarn;
00120
00121 G4double GammaEnergySave = GammaEnergy;
00122 if (GammaEnergy < GammaEnergyLimit) { GammaEnergy = GammaEnergyLimit; }
00123
00124 G4double X=log(GammaEnergy/electron_mass_c2), X2=X*X, X3=X2*X, X4=X3*X, X5=X4*X;
00125
00126 G4double F1 = a0 + a1*X + a2*X2 + a3*X3 + a4*X4 + a5*X5,
00127 F2 = b0 + b1*X + b2*X2 + b3*X3 + b4*X4 + b5*X5,
00128 F3 = c0 + c1*X + c2*X2 + c3*X3 + c4*X4 + c5*X5;
00129
00130 xSection = (Z + 1.)*(F1*Z + F2*Z*Z + F3);
00131
00132 if (GammaEnergySave < GammaEnergyLimit) {
00133
00134 X = (GammaEnergySave - 2.*electron_mass_c2)
00135 / (GammaEnergyLimit - 2.*electron_mass_c2);
00136 xSection *= X*X;
00137 }
00138
00139 if (xSection < 0.) { xSection = 0.; }
00140 return xSection;
00141 }
00142
00143
00144
00145 void G4BetheHeitlerModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
00146 const G4MaterialCutsCouple* couple,
00147 const G4DynamicParticle* aDynamicGamma,
00148 G4double,
00149 G4double)
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162 {
00163 const G4Material* aMaterial = couple->GetMaterial();
00164
00165 G4double GammaEnergy = aDynamicGamma->GetKineticEnergy();
00166 G4ParticleMomentum GammaDirection = aDynamicGamma->GetMomentumDirection();
00167
00168 G4double epsil ;
00169 G4double epsil0 = electron_mass_c2/GammaEnergy ;
00170 if(epsil0 > 1.0) { return; }
00171
00172
00173 static const G4double Egsmall=2.*MeV;
00174
00175
00176 const G4Element* anElement = SelectRandomAtom(aMaterial, theGamma, GammaEnergy);
00177
00178 if (GammaEnergy < Egsmall) {
00179
00180 epsil = epsil0 + (0.5-epsil0)*G4UniformRand();
00181
00182 } else {
00183
00184
00185
00186 G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3());
00187 if (GammaEnergy > 50.*MeV) { FZ += 8.*(anElement->GetfCoulomb()); }
00188
00189
00190 G4double screenfac = 136.*epsil0/(anElement->GetIonisation()->GetZ3());
00191 G4double screenmax = exp ((42.24 - FZ)/8.368) - 0.952 ;
00192 G4double screenmin = min(4.*screenfac,screenmax);
00193
00194
00195 G4double epsil1 = 0.5 - 0.5*sqrt(1. - screenmin/screenmax) ;
00196 G4double epsilmin = max(epsil0,epsil1) , epsilrange = 0.5 - epsilmin;
00197
00198
00199
00200
00201
00202 G4double screenvar, greject ;
00203
00204 G4double F10 = ScreenFunction1(screenmin) - FZ;
00205 G4double F20 = ScreenFunction2(screenmin) - FZ;
00206 G4double NormF1 = max(F10*epsilrange*epsilrange,0.);
00207 G4double NormF2 = max(1.5*F20,0.);
00208
00209 do {
00210 if ( NormF1/(NormF1+NormF2) > G4UniformRand() ) {
00211 epsil = 0.5 - epsilrange*pow(G4UniformRand(), 0.333333);
00212 screenvar = screenfac/(epsil*(1-epsil));
00213 greject = (ScreenFunction1(screenvar) - FZ)/F10;
00214
00215 } else {
00216 epsil = epsilmin + epsilrange*G4UniformRand();
00217 screenvar = screenfac/(epsil*(1-epsil));
00218 greject = (ScreenFunction2(screenvar) - FZ)/F20;
00219 }
00220
00221 } while( greject < G4UniformRand() );
00222
00223 }
00224
00225
00226
00227
00228
00229 G4double ElectTotEnergy, PositTotEnergy;
00230 if (G4UniformRand() > 0.5) {
00231
00232 ElectTotEnergy = (1.-epsil)*GammaEnergy;
00233 PositTotEnergy = epsil*GammaEnergy;
00234
00235 } else {
00236
00237 PositTotEnergy = (1.-epsil)*GammaEnergy;
00238 ElectTotEnergy = epsil*GammaEnergy;
00239 }
00240
00241
00242
00243
00244
00245
00246
00247
00248 G4double u;
00249 const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ;
00250
00251 if (9./(9.+d) >G4UniformRand()) u= - log(G4UniformRand()*G4UniformRand())/a1;
00252 else u= - log(G4UniformRand()*G4UniformRand())/a2;
00253
00254 G4double TetEl = u*electron_mass_c2/ElectTotEnergy;
00255 G4double TetPo = u*electron_mass_c2/PositTotEnergy;
00256 G4double Phi = twopi * G4UniformRand();
00257 G4double dxEl= sin(TetEl)*cos(Phi),dyEl= sin(TetEl)*sin(Phi),dzEl=cos(TetEl);
00258 G4double dxPo=-sin(TetPo)*cos(Phi),dyPo=-sin(TetPo)*sin(Phi),dzPo=cos(TetPo);
00259
00260
00261
00262
00263
00264
00265
00266 G4double ElectKineEnergy = max(0.,ElectTotEnergy - electron_mass_c2);
00267
00268 G4ThreeVector ElectDirection (dxEl, dyEl, dzEl);
00269 ElectDirection.rotateUz(GammaDirection);
00270
00271
00272 G4DynamicParticle* aParticle1= new G4DynamicParticle(
00273 theElectron,ElectDirection,ElectKineEnergy);
00274
00275
00276
00277 G4double PositKineEnergy = max(0.,PositTotEnergy - electron_mass_c2);
00278
00279 G4ThreeVector PositDirection (dxPo, dyPo, dzPo);
00280 PositDirection.rotateUz(GammaDirection);
00281
00282
00283 G4DynamicParticle* aParticle2= new G4DynamicParticle(
00284 thePositron,PositDirection,PositKineEnergy);
00285
00286
00287 fvect->push_back(aParticle1);
00288 fvect->push_back(aParticle2);
00289
00290
00291 fParticleChange->SetProposedKineticEnergy(0.);
00292 fParticleChange->ProposeTrackStatus(fStopAndKill);
00293 }
00294
00295