#include <G4BetheHeitlerModel.hh>
Inheritance diagram for G4BetheHeitlerModel:
Public Member Functions | |
G4BetheHeitlerModel (const G4ParticleDefinition *p=0, const G4String &nam="BetheHeitler") | |
virtual | ~G4BetheHeitlerModel () |
virtual void | Initialise (const G4ParticleDefinition *, const G4DataVector &) |
virtual G4double | ComputeCrossSectionPerAtom (const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0., G4double cut=0., G4double emax=DBL_MAX) |
virtual void | SampleSecondaries (std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy) |
Definition at line 58 of file G4BetheHeitlerModel.hh.
G4BetheHeitlerModel::G4BetheHeitlerModel | ( | const G4ParticleDefinition * | p = 0 , |
|
const G4String & | nam = "BetheHeitler" | |||
) |
Definition at line 68 of file G4BetheHeitlerModel.cc.
References G4Electron::Electron(), G4Gamma::Gamma(), and G4Positron::Positron().
00070 : G4VEmModel(nam) 00071 { 00072 fParticleChange = 0; 00073 theGamma = G4Gamma::Gamma(); 00074 thePositron = G4Positron::Positron(); 00075 theElectron = G4Electron::Electron(); 00076 }
G4BetheHeitlerModel::~G4BetheHeitlerModel | ( | ) | [virtual] |
G4double G4BetheHeitlerModel::ComputeCrossSectionPerAtom | ( | const G4ParticleDefinition * | , | |
G4double | kinEnergy, | |||
G4double | Z, | |||
G4double | A = 0. , |
|||
G4double | cut = 0. , |
|||
G4double | emax = DBL_MAX | |||
) | [virtual] |
Reimplemented from G4VEmModel.
Definition at line 95 of file G4BetheHeitlerModel.cc.
00102 : sigma=sigma(1.5MeV)*(GammaEnergy-2electronmass) 00103 // *(GammaEnergy-2electronmass) 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 }
void G4BetheHeitlerModel::Initialise | ( | const G4ParticleDefinition * | , | |
const G4DataVector & | ||||
) | [virtual] |
Implements G4VEmModel.
Reimplemented in G4PolarizedGammaConversionModel.
Definition at line 85 of file G4BetheHeitlerModel.cc.
References G4VEmModel::GetParticleChangeForGamma(), and G4VEmModel::InitialiseElementSelectors().
Referenced by G4PolarizedGammaConversionModel::Initialise().
00087 { 00088 if(!fParticleChange) { fParticleChange = GetParticleChangeForGamma(); } 00089 InitialiseElementSelectors(p, cuts); 00090 }
void G4BetheHeitlerModel::SampleSecondaries | ( | std::vector< G4DynamicParticle * > * | , | |
const G4MaterialCutsCouple * | , | |||
const G4DynamicParticle * | , | |||
G4double | tmin, | |||
G4double | maxEnergy | |||
) | [virtual] |
Implements G4VEmModel.
Reimplemented in G4PolarizedGammaConversionModel.
Definition at line 145 of file G4BetheHeitlerModel.cc.
References F10, F20, fStopAndKill, G4UniformRand, G4Element::GetfCoulomb(), G4Element::GetIonisation(), G4DynamicParticle::GetKineticEnergy(), G4IonisParamElm::GetlogZ3(), G4MaterialCutsCouple::GetMaterial(), G4DynamicParticle::GetMomentumDirection(), G4IonisParamElm::GetZ3(), G4VParticleChange::ProposeTrackStatus(), G4VEmModel::SelectRandomAtom(), and G4ParticleChangeForGamma::SetProposedKineticEnergy().
Referenced by G4PolarizedGammaConversionModel::SampleSecondaries().
00157 : Effects due to the breakdown of the Born approximation at 00158 // low energy are ignored. 00159 // Note 2 : The differential cross section implicitly takes account of 00160 // pair creation in both nuclear and atomic electron fields. 00161 // However triplet prodution is not generated. 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 // do it fast if GammaEnergy < 2. MeV 00173 static const G4double Egsmall=2.*MeV; 00174 00175 // select randomly one element constituing the material 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 // now comes the case with GammaEnergy >= 2. MeV 00184 00185 // Extract Coulomb factor for this Element 00186 G4double FZ = 8.*(anElement->GetIonisation()->GetlogZ3()); 00187 if (GammaEnergy > 50.*MeV) { FZ += 8.*(anElement->GetfCoulomb()); } 00188 00189 // limits of the screening variable 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 // limits of the energy sampling 00195 G4double epsil1 = 0.5 - 0.5*sqrt(1. - screenmin/screenmax) ; 00196 G4double epsilmin = max(epsil0,epsil1) , epsilrange = 0.5 - epsilmin; 00197 00198 // 00199 // sample the energy rate of the created electron (or positron) 00200 // 00201 //G4double epsil, screenvar, greject ; 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 } // end of epsil sampling 00224 00225 // 00226 // fixe charges randomly 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 // scattered electron (positron) angles. ( Z - axis along the parent photon) 00243 // 00244 // universal distribution suggested by L. Urban 00245 // (Geant3 manual (1993) Phys211), 00246 // derived from Tsai distribution (Rev Mod Phys 49,421(1977)) 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 // kinematic of the created pair 00262 // 00263 // the electron and positron are assumed to have a symetric 00264 // angular distribution with respect to the Z axis along the parent photon. 00265 00266 G4double ElectKineEnergy = max(0.,ElectTotEnergy - electron_mass_c2); 00267 00268 G4ThreeVector ElectDirection (dxEl, dyEl, dzEl); 00269 ElectDirection.rotateUz(GammaDirection); 00270 00271 // create G4DynamicParticle object for the particle1 00272 G4DynamicParticle* aParticle1= new G4DynamicParticle( 00273 theElectron,ElectDirection,ElectKineEnergy); 00274 00275 // the e+ is always created (even with Ekine=0) for further annihilation. 00276 00277 G4double PositKineEnergy = max(0.,PositTotEnergy - electron_mass_c2); 00278 00279 G4ThreeVector PositDirection (dxPo, dyPo, dzPo); 00280 PositDirection.rotateUz(GammaDirection); 00281 00282 // create G4DynamicParticle object for the particle2 00283 G4DynamicParticle* aParticle2= new G4DynamicParticle( 00284 thePositron,PositDirection,PositKineEnergy); 00285 00286 // Fill output vector 00287 fvect->push_back(aParticle1); 00288 fvect->push_back(aParticle2); 00289 00290 // kill incident photon 00291 fParticleChange->SetProposedKineticEnergy(0.); 00292 fParticleChange->ProposeTrackStatus(fStopAndKill); 00293 }