89 lowEnergyLimit (250*
eV),
90 highEnergyLimit(100*
GeV),
91 intrinsicLowEnergyLimit(10*
eV),
92 intrinsicHighEnergyLimit(100*
GeV)
97 G4Exception(
"G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()",
99 "Energy outside intrinsic process validity range!");
106 G4String scatterFile =
"comp/ce-sf-";
143 G4String crossSectionFile =
"comp/ce-cs-";
184 if(!(gammaPolarization0.
isOrthogonal(gammaDirection0, 1
e-6))||(gammaPolarization0.
mag()==0))
221 G4double alpha1 = - std::log(epsilon0);
222 G4double alpha2 = 0.5*(1.- epsilon0Sq);
237 epsilon = std::sqrt(epsilonSq);
240 onecost = (1.-
epsilon)/(epsilon*E0_m);
241 sinThetaSqr = onecost*(2.-onecost);
244 if (sinThetaSqr > 1.)
247 <<
" -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
248 <<
"sin(theta)**2 = "
254 if (sinThetaSqr < 0.)
257 <<
" -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
258 <<
"sin(theta)**2 = "
266 G4double x = std::sqrt(onecost/2.) / (wlGamma/
cm);;
268 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
290 <<
" -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
300 <<
" -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
310 G4double sinTheta = std::sqrt (sinThetaSqr);
316 <<
" -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
326 <<
" -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt "
336 G4double dirx = sinTheta*std::cos(phi);
337 G4double diry = sinTheta*std::sin(phi);
352 G4int maxDopplerIterations = 1000;
354 G4double photonEoriginal = epsilon * gammaEnergy0;
366 eMax = gammaEnergy0 - bindingE;
372 G4double pDoppler2 = pDoppler * pDoppler;
373 G4double var2 = 1. + onecost * E0_m;
374 G4double var3 = var2*var2 - pDoppler2;
375 G4double var4 = var2 - pDoppler2 * cosTheta;
376 G4double var = var4*var4 - var3 + pDoppler2 * var3;
382 if (
G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
383 else photonE = (var4 + varSqrt) * scale;
389 }
while ( iteration <= maxDopplerIterations &&
390 (photonE < 0. || photonE > eMax || photonE < eMax*
G4UniformRand()) );
394 if (iteration >= maxDopplerIterations)
396 photonE = photonEoriginal;
400 gammaEnergy1 = photonE;
426 gammaDirection1 = tmpDirection1;
431 gammaPolarization0,gammaPolarization1);
433 if (gammaEnergy1 > 0.)
449 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
460 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
461 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
496 b = energyRate + 1/energyRate;
498 phiProbability = 1 - (a/
b)*(std::cos(phi)*std::cos(phi));
503 while ( rand2 > phiProbability );
536 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
537 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
538 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
561 return gammaPolarization - gammaPolarization.
dot(gammaDirection)/gammaDirection.
dot(gammaDirection) * gammaDirection;
574 G4double sinTheta = std::sqrt(sinSqrTh);
578 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
617 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
632 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
636 G4double xParallel = normalisation*cosBeta;
637 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
638 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
640 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
641 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
643 G4double xTotal = (xParallel + xPerpendicular);
644 G4double yTotal = (yParallel + yPerpendicular);
645 G4double zTotal = (zParallel + zPerpendicular);
647 gammaPolarization1.
setX(xTotal);
648 gammaPolarization1.
setY(yTotal);
649 gammaPolarization1.
setZ(zTotal);
651 return gammaPolarization1;
672 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
677 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
695 size_t materialIndex = couple->
GetIndex();