125{
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141 if (verboseLevel > 3)
142 G4cout <<
"Calling SamplingSecondaries() of G4PenelopeAnnihilationModel" <<
G4endl;
143
145
146
149
150 if (kineticEnergy == 0.0)
151 {
152
154 G4double sinTheta = std::sqrt(1.0-cosTheta*cosTheta);
156 G4ThreeVector direction (sinTheta*std::cos(phi),sinTheta*std::sin(phi),cosTheta);
158 direction, electron_mass_c2);
160 -direction, electron_mass_c2);
161
162 fvect->push_back(firstGamma);
163 fvect->push_back(secondGamma);
164 return;
165 }
166
167
170 G4double gamma = 1.0 + std::max(kineticEnergy,1.0*eV)/electron_mass_c2;
171 G4double gamma21 = std::sqrt(gamma*gamma-1);
173 G4double chimin = 1.0/(ani+gamma21);
174 G4double rchi = (1.0-chimin)/chimin;
178 do{
180 G4double reject = ani*ani*(1.0-epsilon)+2.0*gamma-(1.0/epsilon);
182 }while(test>0);
183
184 G4double totalAvailableEnergy = kineticEnergy + 2.0*electron_mass_c2;
185 G4double photon1Energy = epsilon*totalAvailableEnergy;
186 G4double photon2Energy = (1.0-epsilon)*totalAvailableEnergy;
187 G4double cosTheta1 = (ani-1.0/epsilon)/gamma21;
188 G4double cosTheta2 = (ani-1.0/(1.0-epsilon))/gamma21;
189
190
191
192 G4double sinTheta1 = std::sqrt(1.-cosTheta1*cosTheta1);
194 G4double dirx1 = sinTheta1 * std::cos(phi1);
195 G4double diry1 = sinTheta1 * std::sin(phi1);
197
198 G4double sinTheta2 = std::sqrt(1.-cosTheta2*cosTheta2);
200 G4double dirx2 = sinTheta2 * std::cos(phi2);
201 G4double diry2 = sinTheta2 * std::sin(phi2);
203
205 photon1Direction.rotateUz(positronDirection);
206
208 photon1Direction,
209 photon1Energy);
210 fvect->push_back(aParticle1);
211
213 photon2Direction.rotateUz(positronDirection);
214
216 photon2Direction,
217 photon2Energy);
218 fvect->push_back(aParticle2);
219
220 if (verboseLevel > 1)
221 {
222 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
223 G4cout <<
"Energy balance from G4PenelopeAnnihilation" <<
G4endl;
224 G4cout <<
"Kinetic positron energy: " << kineticEnergy/keV <<
" keV" <<
G4endl;
225 G4cout <<
"Total available energy: " << totalAvailableEnergy/keV <<
" keV " <<
G4endl;
226 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
227 G4cout <<
"Photon energy 1: " << photon1Energy/keV <<
" keV" <<
G4endl;
228 G4cout <<
"Photon energy 2: " << photon2Energy/keV <<
" keV" <<
G4endl;
229 G4cout <<
"Total final state: " << (photon1Energy+photon2Energy)/keV <<
231 G4cout <<
"-----------------------------------------------------------" <<
G4endl;
232 }
233 if (verboseLevel > 0)
234 {
235 G4double energyDiff = std::fabs(totalAvailableEnergy-photon1Energy-photon2Energy);
236 if (energyDiff > 0.05*keV)
237 G4cout <<
"Warning from G4PenelopeAnnihilation: problem with energy conservation: " <<
238 (photon1Energy+photon2Energy)/keV <<
239 " keV (final) vs. " <<
240 totalAvailableEnergy/keV <<
" keV (initial)" <<
G4endl;
241 }
242 return;
243}
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void ProposeTrackStatus(G4TrackStatus status)