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
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069
00070
00071 #include "G4eeToTwoGammaModel.hh"
00072 #include "G4PhysicalConstants.hh"
00073 #include "G4SystemOfUnits.hh"
00074 #include "G4TrackStatus.hh"
00075 #include "G4Electron.hh"
00076 #include "G4Positron.hh"
00077 #include "G4Gamma.hh"
00078 #include "Randomize.hh"
00079 #include "G4ParticleChangeForGamma.hh"
00080
00081
00082
00083 using namespace std;
00084
00085 G4eeToTwoGammaModel::G4eeToTwoGammaModel(const G4ParticleDefinition*,
00086 const G4String& nam)
00087 : G4VEmModel(nam),
00088 pi_rcl2(pi*classic_electr_radius*classic_electr_radius),
00089 isInitialised(false)
00090 {
00091 theGamma = G4Gamma::Gamma();
00092 fParticleChange = 0;
00093 }
00094
00095
00096
00097 G4eeToTwoGammaModel::~G4eeToTwoGammaModel()
00098 {}
00099
00100
00101
00102 void G4eeToTwoGammaModel::Initialise(const G4ParticleDefinition*,
00103 const G4DataVector&)
00104 {
00105 if(isInitialised) { return; }
00106 fParticleChange = GetParticleChangeForGamma();
00107 isInitialised = true;
00108 }
00109
00110
00111
00112 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerElectron(
00113 const G4ParticleDefinition*,
00114 G4double kineticEnergy,
00115 G4double, G4double)
00116 {
00117
00118
00119
00120 G4double ekin = std::max(eV,kineticEnergy);
00121
00122 G4double tau = ekin/electron_mass_c2;
00123 G4double gam = tau + 1.0;
00124 G4double gamma2= gam*gam;
00125 G4double bg2 = tau * (tau+2.0);
00126 G4double bg = sqrt(bg2);
00127
00128 G4double cross = pi_rcl2*((gamma2+4*gam+1.)*log(gam+bg) - (gam+3.)*bg)
00129 / (bg2*(gam+1.));
00130 return cross;
00131 }
00132
00133
00134
00135 G4double G4eeToTwoGammaModel::ComputeCrossSectionPerAtom(
00136 const G4ParticleDefinition* p,
00137 G4double kineticEnergy, G4double Z,
00138 G4double, G4double, G4double)
00139 {
00140
00141
00142 G4double cross = Z*ComputeCrossSectionPerElectron(p,kineticEnergy);
00143 return cross;
00144 }
00145
00146
00147
00148 G4double G4eeToTwoGammaModel::CrossSectionPerVolume(
00149 const G4Material* material,
00150 const G4ParticleDefinition* p,
00151 G4double kineticEnergy,
00152 G4double, G4double)
00153 {
00154
00155
00156 G4double eDensity = material->GetElectronDensity();
00157 G4double cross = eDensity*ComputeCrossSectionPerElectron(p,kineticEnergy);
00158 return cross;
00159 }
00160
00161
00162
00163 void G4eeToTwoGammaModel::SampleSecondaries(vector<G4DynamicParticle*>* vdp,
00164 const G4MaterialCutsCouple*,
00165 const G4DynamicParticle* dp,
00166 G4double,
00167 G4double)
00168 {
00169 G4double PositKinEnergy = dp->GetKineticEnergy();
00170 G4DynamicParticle *aGamma1, *aGamma2;
00171
00172
00173 if(PositKinEnergy == 0.0) {
00174 G4double cost = 2.*G4UniformRand()-1.;
00175 G4double sint = sqrt((1. - cost)*(1. + cost));
00176 G4double phi = twopi * G4UniformRand();
00177 G4ThreeVector dir(sint*cos(phi), sint*sin(phi), cost);
00178 phi = twopi * G4UniformRand();
00179 G4ThreeVector pol(cos(phi), sin(phi), 0.0);
00180 pol.rotateUz(dir);
00181 aGamma1 = new G4DynamicParticle(theGamma, dir, electron_mass_c2);
00182 aGamma1->SetPolarization(pol.x(),pol.y(),pol.z());
00183 aGamma2 = new G4DynamicParticle(theGamma,-dir, electron_mass_c2);
00184 aGamma1->SetPolarization(-pol.x(),-pol.y(),-pol.z());
00185
00186 } else {
00187
00188 G4ThreeVector PositDirection = dp->GetMomentumDirection();
00189
00190 G4double tau = PositKinEnergy/electron_mass_c2;
00191 G4double gam = tau + 1.0;
00192 G4double tau2 = tau + 2.0;
00193 G4double sqgrate = sqrt(tau/tau2)*0.5;
00194 G4double sqg2m1 = sqrt(tau*tau2);
00195
00196
00197 G4double epsilmin = 0.5 - sqgrate;
00198 G4double epsilmax = 0.5 + sqgrate;
00199 G4double epsilqot = epsilmax/epsilmin;
00200
00201
00202
00203
00204 G4double epsil, greject;
00205
00206 do {
00207 epsil = epsilmin*pow(epsilqot,G4UniformRand());
00208 greject = 1. - epsil + (2.*gam*epsil-1.)/(epsil*tau2*tau2);
00209 } while( greject < G4UniformRand() );
00210
00211
00212
00213
00214
00215 G4double cost = (epsil*tau2-1.)/(epsil*sqg2m1);
00216 if(std::abs(cost) > 1.0) {
00217 G4cout << "### G4eeToTwoGammaModel WARNING cost= " << cost
00218 << " positron Ekin(MeV)= " << PositKinEnergy
00219 << " gamma epsil= " << epsil
00220 << G4endl;
00221 if(cost > 1.0) cost = 1.0;
00222 else cost = -1.0;
00223 }
00224 G4double sint = sqrt((1.+cost)*(1.-cost));
00225 G4double phi = twopi * G4UniformRand();
00226
00227
00228
00229
00230
00231 G4double TotalAvailableEnergy = PositKinEnergy + 2.0*electron_mass_c2;
00232 G4double Phot1Energy = epsil*TotalAvailableEnergy;
00233
00234 G4ThreeVector Phot1Direction(sint*cos(phi), sint*sin(phi), cost);
00235 Phot1Direction.rotateUz(PositDirection);
00236 aGamma1 = new G4DynamicParticle (theGamma,Phot1Direction, Phot1Energy);
00237 phi = twopi * G4UniformRand();
00238 G4ThreeVector pol1(cos(phi), sin(phi), 0.0);
00239 pol1.rotateUz(Phot1Direction);
00240 aGamma1->SetPolarization(pol1.x(),pol1.y(),pol1.z());
00241
00242 G4double Phot2Energy =(1.-epsil)*TotalAvailableEnergy;
00243 G4double PositP= sqrt(PositKinEnergy*(PositKinEnergy+2.*electron_mass_c2));
00244 G4ThreeVector dir = PositDirection*PositP - Phot1Direction*Phot1Energy;
00245 G4ThreeVector Phot2Direction = dir.unit();
00246
00247
00248 aGamma2 = new G4DynamicParticle (theGamma,Phot2Direction, Phot2Energy);
00249
00251 aGamma2->SetPolarization(-pol1.x(),-pol1.y(),-pol1.z());
00252 }
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262 vdp->push_back(aGamma1);
00263 vdp->push_back(aGamma2);
00264 fParticleChange->SetProposedKineticEnergy(0.);
00265 fParticleChange->ProposeTrackStatus(fStopAndKill);
00266 }
00267
00268