74 const G4int arrayDim = 980;
84 for(
G4int level = 0; level < 2; level++){
86 char nameChar0[100] =
"ftab0.dat";
87 char nameChar1[100] =
"ftab1.dat";
90 if(level == 0) filename = nameChar0;
91 if(level == 1) filename = nameChar1;
93 char* path = getenv(
"G4LEDATA");
96 G4String excep =
"G4EMDataSet - G4LEDATA environment variable not set";
97 G4Exception(
"G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
103 G4String dirFile = pathString +
"/photoelectric_angular/" + filename;
104 std::ifstream infile(dirFile);
105 if (!infile.is_open())
107 G4String excep =
"data file: " + dirFile +
" not found";
108 G4Exception(
"G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized",
115 for(
G4int i=0 ; i<arrayDim ;i++){
117 infile >>
beta >> aRead >> cRead;
140 G4double beta = std::sqrt((gamma - 1)*(gamma + 1))/gamma;
178 G4double crossSectionMajorantFunctionValue = 0;
195 theta=std::sqrt(((
G4Exp(rand2*std::log(1+cBeta*
pi*
pi)))-1)/cBeta);
196 crossSectionMajorantFunctionValue =
203 theta = std::sqrt(((
G4Exp(rand2*std::log(1+cBeta*
pi*
pi)))-1)/cBeta);
204 crossSectionMajorantFunctionValue =
210 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
213 if(crossSectionValue < 0.0) { crossSectionValue = maxBeta; }
215 }
while(maxBeta > crossSectionValue || theta >
CLHEP::pi);
226 G4double crossSectionMajorantFunctionValue = 0;
227 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
228 return crossSectionMajorantFunctionValue;
239 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
240 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
241 G4double cosTheta = std::cos(theta);
242 G4double sinTheta2 = std::sin(theta)*std::sin(theta);
243 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
244 G4double oneBetaCosTheta = 1-beta*cosTheta;
249 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
250 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
251 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
253 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
254 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
255 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
256 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
257 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
258 (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
275 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
276 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
277 G4double cosTheta = std::cos(theta);
278 G4double sinTheta2 =std::sin(theta)*std::sin(theta);
279 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
280 G4double oneBetaCosTheta = 1-beta*cosTheta;
286 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
287 * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
288 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
290 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
291 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
292 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
293 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
294 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
295 (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
312 if(!(polarization.
isOrthogonal(direction,kTolerance)) || mS == 0){
320 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
321 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
322 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
323 polarization2 = c.
unit();
324 mS = polarization2.
mag();
329 polarization2 = polarization
330 - polarization.
dot(direction)/direction.
dot(direction) * direction;
335 polarization2 = polarization2/mS;
353 if(shellId > 0) { level = 1; }
378 else if(k == indexMax)
385 *majorantSurfaceParameterA = aBeta;
386 *majorantSurfaceParameterC = cBeta;
399 G4ThreeVector outgoingDirection = rotation*samplingDirection;
400 return outgoingDirection;
406 G4cout <<
"Polarized Photoelectric Angular Generator" <<
G4endl;
407 G4cout <<
"PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" <<
G4endl;
408 G4cout <<
"Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." <<
G4endl;
409 G4cout <<
"For higher shells the L1 cross-section is used." <<
G4endl;
410 G4cout <<
"(see Physics Reference Manual) \n" <<
G4endl;
G4ThreeVector PhotoElectronComputeFinalDirection(const G4RotationMatrix &rotation, G4double theta, G4double phi) const
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void PhotoElectronGetMajorantSurfaceAandCParameters(G4int shellId, G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
G4RotationMatrix PhotoElectronRotationMatrix(const G4ThreeVector &direction, const G4ThreeVector &polarization)
G4double CrossSectionMajorantFunction(G4double theta, G4double cBeta) const
const G4ThreeVector & GetMomentumDirection() const
G4double aMajorantSurfaceParameterTable[980][2]
double dot(const Hep3Vector &) const
G4double cMajorantSurfaceParameterTable[980][2]
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
G4ThreeVector PerpendicularVector(const G4ThreeVector &a) const
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double eKinEnergy, G4int shellId, const G4Material *mat=0)
static constexpr double electron_mass_c2
G4double DSigmaL1shellGavrila(G4double beta, G4double theta, G4double phi) const
static constexpr double twopi
G4double DSigmaKshellGavrila1959(G4double beta, G4double theta, G4double phi) const
double howOrthogonal(const Hep3Vector &v) const
Hep3Vector cross(const Hep3Vector &) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4PhotoElectricAngularGeneratorPolarized()
const G4ThreeVector & GetPolarization() const
G4ThreeVector fLocalDirection
G4GLOB_DLL std::ostream G4cout
void PhotoElectronGeneratePhiAndTheta(G4int shellId, G4double beta, G4double aBeta, G4double cBeta, G4double *pphi, G4double *ptheta) const
static constexpr double pi
static constexpr double fine_structure_const
~G4PhotoElectricAngularGeneratorPolarized()
void PrintGeneratorInformation() const
static constexpr double pi