272{
274
275
279 tmax = 0.5*kineticEnergy;
280 } else {
281 tmax = kineticEnergy;
282 }
283 if(maxEnergy < tmax) { tmax = maxEnergy; }
284 if(tmin >= tmax) { return; }
285
295
296
298
301 grej = 1.0 - gg*xmax + xmax*xmax*(1.0 - gg + (1.0 - gg*y)/(y*y));
302
303 do {
305 x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
306 y = 1.0 - x;
307 z = 1.0 - gg*x + x*x*(1.0 - gg + (1.0 - gg*y)/(y*y));
308
309
310
311
312
313
314
315
316
317
318 } while(grej * rndm[1] > z);
319
320
321 } else {
322
331
332 y = xmax*xmax;
333 grej = 1.0 + (y*y*b4 - xmin*xmin*xmin*b3 + y*b2 - xmin*b1)*beta2;
334 do {
336 x = xmin*xmax/(xmin*(1.0 - rndm[0]) + xmax*rndm[0]);
337 y = x*x;
338 z = 1.0 + (y*y*b4 - x*y*b3 + y*b2 - x*b1)*beta2;
339
340
341
342
343
344
345
346
347
348
349 } while(grej * rndm[1] > z);
350 }
351
352 G4double deltaKinEnergy = x * kineticEnergy;
353
355
359
360 deltaDirection =
362
363 } else {
364
366 sqrt(deltaKinEnergy * (deltaKinEnergy + 2.0*electron_mass_c2));
369 if(cost > 1.0) { cost = 1.0; }
370 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
371
373
374 deltaDirection.
set(sint*cos(phi),sint*sin(phi), cost) ;
376 }
377
378
380 vdp->push_back(delta);
381
382
383 kineticEnergy -= deltaKinEnergy;
385 finalP = finalP.
unit();
386
389}
void set(double x, double y, double z)
Hep3Vector & rotateUz(const Hep3Vector &)
virtual void flatArray(const int size, double *vect)=0
const G4ThreeVector & GetMomentumDirection() const
G4double GetKineticEnergy() const
G4ThreeVector GetMomentum() const
G4double GetTotalMomentum() const
const G4Material * GetMaterial() const
void SetProposedKineticEnergy(G4double proposedKinEnergy)
void SetProposedMomentumDirection(const G4ThreeVector &dir)
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
G4int SelectRandomAtomNumber(const G4Material *) const
G4double energy(const ThreeVector &p, const G4double m)