157{
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173 if (verboseLevel > 3)
174 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" <<
G4endl;
175
177
178
181
182 if (kineticEnergy == 0.0)
183 {
184
186 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
188 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
190 direction, electron_mass_c2);
192 -direction, electron_mass_c2);
193
194 fvect->push_back(firstGamma);
195 fvect->push_back(secondGamma);
196 return;
197 }
198
199
202 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
203 G4double gamma21 = std::sqrt(gamma*gamma-1);
205 G4double chimin = 1.0/(ani+gamma21);
206 G4double rchi = (1.0-chimin)/chimin;
210 do{
214 }while(test>0);
215
216 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
221
222
223
224 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
226 G4double dirx1 = sinTheta1 * std::cos(phi1);
227 G4double diry1 = sinTheta1 * std::sin(phi1);
229
230 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
232 G4double dirx2 = sinTheta2 * std::cos(phi2);
233 G4double diry2 = sinTheta2 * std::sin(phi2);
235
237 photon1Direction.rotateUz(positronDirection);
238
240 photon1Direction,
241 photon1Energy);
242 fvect->push_back(aParticle1);
243
245 photon2Direction.rotateUz(positronDirection);
246
248 photon2Direction,
249 photon2Energy);
250 fvect->push_back(aParticle2);
251
252 if (verboseLevel > 1)
253 {
254 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
255 G4cout <<
"Energy balance from G4PenelopeAnnihilation" <<
G4endl;
256 G4cout <<
"Kinetic positron energy: " << kineticEnergy/keV <<
" keV" <<
G4endl;
257 G4cout <<
"Total available energy: " << totalAvailableEnergy/keV <<
" keV " <<
G4endl;
258 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
259 G4cout <<
"Photon energy 1: " << photon1Energy/keV <<
" keV" <<
G4endl;
260 G4cout <<
"Photon energy 2: " << photon2Energy/keV <<
" keV" <<
G4endl;
261 G4cout <<
"Total final state: " << (photon1Energy+photon2Energy)/keV <<
263 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
264 }
265 if (verboseLevel > 0)
266 {
267 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
268 if (energyDiff > 0.05*keV)
269 G4cout <<
"Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
270 (photon1Energy+photon2Energy)/keV <<
271 " keV (final) vs. " <<
272 totalAvailableEnergy/keV <<
" keV (initial)" <<
G4endl;
273 }
274 return;
275}
double epsilon(double density, double temperature)
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeTrackStatus(G4TrackStatus status)