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
281 do {
282 cost =
G4Exp(logcostm*rndmEngine->
flat());
284 G4double bre = (1.-5.*cost*cost)/(2.*cost);
285 loga =
G4Log((1.+ cost)/(1.- cost));
286 f1_re = 1. - bre*loga;
287 greject = (cost < costlim) ? are*f1_re : 1.0;
288 }
while(greject < rndmEngine->
flat());
289
290
291 G4double sint2 = (1. - cost)*(1. + cost);
292 G4double fp = 1. - sint2*loga/(2.*cost) ;
294 do {
295 phi_re = twopi*rndmEngine->
flat();
296 rt = (1. - std::cos(2.*phi_re)*fp/f1_re)/twopi;
297 }
while(rt < rndmEngine->
flat());
298
299
300 G4double S = electron_mass_c2*(2.* photonEnergy + electron_mass_c2);
301 G4double P2 =
S - electron_mass_c2*electron_mass_c2;
302
303 G4double D2 = 4.*
S * electron_mass_c2*electron_mass_c2 + P2*P2*sint2;
304 G4double ener_re = electron_mass_c2 * (
S + electron_mass_c2*electron_mass_c2)/sqrt(D2);
305
306 if(ener_re >= energyThreshold)
307 {
308 G4double electronRKineEnergy = ener_re - electron_mass_c2;
310 G4ThreeVector electronRDirection (sint*std::cos(phi_re), sint*std::sin(phi_re), cost);
311 electronRDirection.rotateUz(photonDirection);
313 electronRDirection,
314 electronRKineEnergy);
315 }
316 else
317 {
318
319
321 ener_re = 0.0;
322 }
323
324
325
326
328
329 G4double a = std::sqrt(16./xb - 3. - 36.*re*xn + 36.*re*re*xn*xn + 6.*xb*re*xn);
331 epsilon = c1/(2.*xb) + (xb - 4.)/(2.*c1) + 0.5;
332
333 G4double photonEnergy1 = photonEnergy - ener_re ;
334
335 G4double positronTotEnergy = std::max(
epsilon*photonEnergy1, electron_mass_c2);
336 G4double electronTotEnergy = std::max(photonEnergy1 - positronTotEnergy, electron_mass_c2);
337
339 static const G4double a2 = 0.5333333333;
341 G4double u = (0.25 > rndmEngine->
flat()) ? uu*a1 : uu*a2;
342
343 G4double thetaEle = u*electron_mass_c2/electronTotEnergy;
344 G4double sinte = std::sin(thetaEle);
345 G4double coste = std::cos(thetaEle);
346
347 G4double thetaPos = u*electron_mass_c2/positronTotEnergy;
348 G4double sintp = std::sin(thetaPos);
349 G4double costp = std::cos(thetaPos);
350
354
355
356
357
358
359 G4double electronKineEnergy = electronTotEnergy - electron_mass_c2;
360
361 G4ThreeVector electronDirection (sinte*cosp, sinte*sinp, coste);
362 electronDirection.rotateUz(photonDirection);
363
365 electronDirection,
366 electronKineEnergy);
367
368 G4double positronKineEnergy = positronTotEnergy - electron_mass_c2;
369
370 G4ThreeVector positronDirection (-sintp*cosp, -sintp*sinp, costp);
371 positronDirection.rotateUz(photonDirection);
372
373
375 positronDirection, positronKineEnergy);
376
377
378 fvect->push_back(particle1);
379 fvect->push_back(particle2);
380
381 if(particle3) { fvect->push_back(particle3); }
382
383
386}
G4double epsilon(G4double density, G4double temperature)
G4double S(G4double temp)
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)