148{
149
150
151
152
153 if (verboseLevel > 3)
154 G4cout <<
"Calling SampleSecondaries() of G4BoldyshevTripletModel" <<
G4endl;
155
158
161
162 G4double positronTotEnergy, electronTotEnergy, thetaEle, thetaPos;
163 G4double ener_re=0., theta_re, phi_re, phi;
164
165
166
167 G4double energyThreshold = sqrt(2.)*electron_mass_c2;
168 energyThreshold = 1.1*electron_mass_c2;
169
170
171 G4double momentumThreshold_c = sqrt(energyThreshold * energyThreshold - electron_mass_c2*electron_mass_c2);
172 G4double momentumThreshold_N = momentumThreshold_c/electron_mass_c2;
173
174
175
176 G4double SigmaTot = (28./9.) * std::log ( 2.* photonEnergy / electron_mass_c2 ) - 218. / 27. ;
177 G4double X_0 = 2. * ( sqrt(momentumThreshold_N*momentumThreshold_N + 1) -1 );
178 G4double SigmaQ = (82./27. - (14./9.) * log (X_0) + 4./15.*X_0 - 0.0348 * X_0 * X_0);
180
181
182 if (recoilProb >= SigmaQ/SigmaTot)
183 {
184
185 G4double cosThetaMax = ( ( energyThreshold - electron_mass_c2 ) / (momentumThreshold_c) + electron_mass_c2*
186 ( energyThreshold + electron_mass_c2 ) / (photonEnergy*momentumThreshold_c) );
187
189
192 G4double are, bre, loga, f1_re, greject, cost;
193
194 do {
197
198 cost = pow(cosThetaMax,r1);
199 theta_re = acos(cost);
200 are = 1./(14.*cost*cost);
201 bre = (1.-5.*cost*cost)/(2.*cost);
202 loga = log((1.+ cost)/(1.- cost));
203 f1_re = 1. - bre*loga;
204
205 if ( theta_re >= 4.47*CLHEP::pi/180.)
206 {
207 greject = are*f1_re;
208 } else {
209 greject = 1. ;
210 }
211 } while(greject < r2);
212
213
214
216
217 do {
218
221 phi_re = twopi*r3 ;
223 G4double fp = 1. - sint2*loga/(2.*cost) ;
224 rt = (1.-cos(2.*phi_re)*fp/f1_re)/(2.*pi) ;
225
226 } while(rt < r4);
227
228
229
230 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
231 G4double D2 = 4.*S * electron_mass_c2*electron_mass_c2
232 + (S - electron_mass_c2*electron_mass_c2)
233 *(S - electron_mass_c2*electron_mass_c2)*sin(theta_re)*sin(theta_re);
234 ener_re = electron_mass_c2 * (S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
235
236
237
238
239 G4double dxEle_re=sin(theta_re)*std::cos(phi_re),dyEle_re=sin(theta_re)*std::sin(phi_re), dzEle_re=cos(theta_re);
240
241 G4double electronRKineEnergy = std::max(0.,ener_re - electron_mass_c2) ;
242
243 G4ThreeVector electronRDirection (dxEle_re, dyEle_re, dzEle_re);
244 electronRDirection.rotateUz(photonDirection);
245
247 electronRDirection,
248 electronRKineEnergy);
249 fvect->push_back(particle3);
250
251 }
252 else
253 {
254
255
257 }
258
259
260
261
262 G4double t = 0.5*log(momentumThreshold_N + sqrt(momentumThreshold_N*momentumThreshold_N+1));
263
264 G4cout << 0.5*asinh(momentumThreshold_N) <<
" " << t <<
G4endl;
265
266 G4double J1 = 0.5*(t*cosh(t)/sinh(t) - log(2.*sinh(t)));
267 G4double J2 = (-2./3.)*log(2.*sinh(t)) + t*cosh(t)/sinh(t) + (sinh(t)-t*pow(cosh(t),3))/(3.*pow(sinh(t),3));
269
274 G4double b1 = 16. - 3.*b - 36.*b*re*
n + 36.*b*pow(re,2.)*pow(n,2.) +
276 a = pow((b1/b),0.5);
277 G4double c1 = (-6. + 12.*re*
n + b + 2*a)*pow(b,2.);
278 epsilon = (pow(c1,1./3.))/(2.*b) + (b-4.)/(2.*pow(c1,1./3.))+0.5;
279
280 G4double photonEnergy1 = photonEnergy - ener_re ;
281 positronTotEnergy = epsilon*photonEnergy1;
282 electronTotEnergy = photonEnergy1 - positronTotEnergy;
283
284 G4double momento_e = sqrt(electronTotEnergy*electronTotEnergy -
285 electron_mass_c2*electron_mass_c2) ;
286 G4double momento_p = sqrt(positronTotEnergy*positronTotEnergy -
287 electron_mass_c2*electron_mass_c2) ;
288
289 thetaEle = acos((sqrt(p0*p0/(momento_e*momento_e) +1.)- p0/momento_e)) ;
290 thetaPos = acos((sqrt(p0*p0/(momento_p*momento_p) +1.)- p0/momento_p)) ;
292
293 G4double dxEle= std::sin(thetaEle)*std::cos(phi),dyEle= std::sin(thetaEle)*std::sin(phi),dzEle=std::cos(thetaEle);
294 G4double dxPos=-std::sin(thetaPos)*std::cos(phi),dyPos=-std::sin(thetaPos)*std::sin(phi),dzPos=std::cos(thetaPos);
295
296
297
298
299
300
301 G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
302
303
304
306 electronDirection.rotateUz(photonDirection);
307
309 electronDirection,
310 electronKineEnergy);
311
312
313 G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
314
315
316
318 positronDirection.rotateUz(photonDirection);
319
320
322 positronDirection, positronKineEnergy);
323
324
325
326 fvect->push_back(particle1);
327 fvect->push_back(particle2);
328
329
330
331
332
335
336}
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
static G4Electron * Electron()
void SetProposedKineticEnergy(G4double proposedKinEnergy)
static G4Positron * Positron()
void ProposeTrackStatus(G4TrackStatus status)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)