250{
251
252
253
254 if (verboseLevel > 1) {
255 G4cout <<
"Calling SampleSecondaries() of G4BoldyshevTripletModel"
257 }
258
261
263
265
266
268 static const G4double costlim = std::cos(4.47*CLHEP::pi/180.);
269
270 G4double loga, f1_re, greject, cost;
271 G4double cosThetaMax = (energyThreshold - electron_mass_c2
272 + electron_mass_c2*(energyThreshold + electron_mass_c2)/photonEnergy )
273 /momentumThreshold_c;
274 if (cosThetaMax > 1.) {
275
276
277 cosThetaMax = 1.0;
278 }
279
282 do {
283 cost =
G4Exp(logcostm*rndmEngine->
flat());
285 G4double bre = (1.-5.*cost*cost)/(2.*cost);
286 loga =
G4Log((1.+ cost)/(1.- cost));
287 f1_re = 1. - bre*loga;
288 greject = (cost < costlim) ? are*f1_re : 1.0;
289
291 }
while(greject < rndmEngine->
flat());
292
293
294 G4double sint2 = (1. - cost)*(1. + cost);
295 G4double fp = 1. - sint2*loga/(2.*cost) ;
298 do {
299 phi_re = twopi*rndmEngine->
flat();
300 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi;
301
303 }
while(rt < rndmEngine->
flat());
304
305
306 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
307 G4double P2 =
S - electron_mass_c2*electron_mass_c2;
308
309 G4double D2 = 4.*
S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2;
310 G4double ener_re = electron_mass_c2 * (
S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
311
312 if(ener_re >= energyThreshold)
313 {
314 G4double electronRKineEnergy = ener_re - electron_mass_c2;
316 G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost);
317 electronRDirection.rotateUz(photonDirection);
319 electronRDirection,
320 electronRKineEnergy);
321 }
322 else
323 {
324
325
327 ener_re = 0.0;
328 }
329
330
331
332
334
335 G4double a = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn);
337 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5;
338
339 G4double photonEnergy1 = photonEnergy - ener_re ;
340
341 G4double positronTotEnergy = std::max(
epsilon*photonEnergy1, electron_mass_c2);
342 G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2);
343
345 static const G4double a2 = 0.5333333333;
347 G4double u = (0.25 > rndmEngine->
flat()) ? uu*a1 : uu*a2;
348
349 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
350 G4double sinte = std::sin(thetaEle);
351 G4double coste = std::cos(thetaEle);
352
353 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
354 G4double sintp = std::sin(thetaPos);
355 G4double costp = std::cos(thetaPos);
356
360
361
362
363
364
365 G4double electronKineEnergy = electronTotEnergy - electron_mass_c2;
366
367 G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
368 electronDirection.rotateUz(photonDirection);
369
371 electronDirection,
372 electronKineEnergy);
373
374 G4double positronKineEnergy = positronTotEnergy - electron_mass_c2;
375
376 G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
377 positronDirection.rotateUz(photonDirection);
378
379
381 positronDirection, positronKineEnergy);
382
383
384 fvect->push_back(particle1);
385 fvect->push_back(particle2);
386
387 if(particle3) { fvect->push_back(particle3); }
388
389
392}
double epsilon(double density, double temperature)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
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)