73 const G4int arrayDim = 980;
83 for(
G4int level = 0; level < 2; level++){
85 char nameChar0[100] =
"ftab0.dat";
86 char nameChar1[100] =
"ftab1.dat";
89 if(level == 0) filename = nameChar0;
90 if(level == 1) filename = nameChar1;
92 char* path = getenv(
"G4LEDATA");
95 G4String excep =
"G4LEDATA environment variable not set!";
96 G4Exception(
"G4RDPhotoElectricAngularGeneratorPolarized()",
101 G4String dirFile = pathString +
"/photoelectric_angular/" + filename;
103 infile = fopen(dirFile,
"r");
106 G4String excep =
"Data file: " + dirFile +
" not found";
107 G4Exception(
"G4RDPhotoElectricAngularGeneratorPolarized()",
113 for(
G4int i=0 ; i<arrayDim ;i++){
114 fscanf(infile,
"%f\t %e\t %e",&beta,&aRead,&cRead);
141 G4int shellLevel = 0;
142 if(shellId < 2) shellLevel = 0;
143 if(shellId >= 2) shellLevel = 1;
157 return final_direction;
170 G4double crossSectionMajorantFunctionValue = 0;
184 theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*
pi*
pi)))-1)/cBeta);
191 theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*
pi*
pi)))-1)/cBeta);
197 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
199 }
while(maxBeta > crossSectionValue);
210 G4double crossSectionMajorantFunctionValue = 0;
211 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
212 return crossSectionMajorantFunctionValue;
224 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
225 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
226 G4double cosTheta = std::cos(theta);
227 G4double sinTheta2 = std::sin(theta)*std::sin(theta);
228 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
229 G4double oneBetaCosTheta = 1-beta*cosTheta;
234 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
235 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
236 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
238 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
239 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
240 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
241 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
242 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
243 (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
259 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
260 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
261 G4double cosTheta = std::cos(theta);
262 G4double sinTheta2 =std::sin(theta)*std::sin(theta);
263 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
264 G4double oneBetaCosTheta = 1-beta*cosTheta;
270 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
271 * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
272 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
274 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
275 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
276 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
277 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
278 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
279 (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
304 if(!(polarization.
isOrthogonal(direction,kTolerance)) || mS == 0){
312 c.
setX(std::cos(angle)*(a0.
x())+std::sin(angle)*b0.
x());
313 c.
setY(std::cos(angle)*(a0.
y())+std::sin(angle)*b0.
y());
314 c.
setZ(std::cos(angle)*(a0.
z())+std::sin(angle)*b0.
z());
315 polarization2 = c.
unit();
316 mS = polarization2.
mag();
321 polarization2 = polarization - polarization.
dot(direction)/direction.
dot(direction) * direction;
326 polarization2 = polarization2/mS;
341 G4int level = shellLevel;
342 if(shellLevel > 1) level = 1;
367 else if(k == indexMax)
374 *majorantSurfaceParameterA = aBeta;
375 *majorantSurfaceParameterC = cBeta;
385 G4double px = std::cos(phi)*std::sin(theta);
386 G4double py = std::sin(phi)*std::sin(theta);
391 G4ThreeVector outgoingDirection = rotation*samplingDirection;
392 return outgoingDirection;
400 G4cout <<
"Polarized Photoelectric Angular Generator" <<
G4endl;
401 G4cout <<
"PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" <<
G4endl;
402 G4cout <<
"Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." <<
G4endl;
403 G4cout <<
"For higher shells the L1 cross-section is used." <<
G4endl;
404 G4cout <<
"(see Physics Reference Manual) \n" <<
G4endl;
void PhotoElectronGeneratePhiAndTheta(const G4int shellLevel, const G4double beta, const G4double aBeta, const G4double cBeta, G4double *pphi, G4double *ptheta) const
CLHEP::Hep3Vector G4ThreeVector
std::vector< ExP01TrackerHit * > a
~G4RDPhotoElectricAngularGeneratorPolarized()
double dot(const Hep3Vector &) const
G4double DSigmaKshellGavrila1959(const G4double beta, const G4double theta, const G4double phi) const
G4double CrossSectionMajorantFunction(const G4double theta, const G4double cBeta) const
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
void PrintGeneratorInformation() const
static constexpr double electron_mass_c2
static constexpr double twopi
G4ThreeVector GetPhotoElectronDirection(const G4ThreeVector &direction, const G4double kineticEnergy, const G4ThreeVector &polarization, const G4int shellId) const
double howOrthogonal(const Hep3Vector &v) const
G4double GetMax(const G4double arg1, const G4double arg2) const
G4double DSigmaL1shellGavrila(const G4double beta, const G4double theta, const G4double phi) const
Hep3Vector cross(const Hep3Vector &) const
G4RDPhotoElectricAngularGeneratorPolarized(const G4String &name)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
void PhotoElectronGetMajorantSurfaceAandCParameters(const G4int shellLevel, const G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
G4RotationMatrix PhotoElectronRotationMatrix(const G4ThreeVector &direction, const G4ThreeVector &polarization) const
G4ThreeVector SetPerpendicularVector(const G4ThreeVector &a) const
G4GLOB_DLL std::ostream G4cout
static constexpr double pi
static constexpr double fine_structure_const
G4ThreeVector PhotoElectronComputeFinalDirection(const G4RotationMatrix &rotation, const G4double theta, const G4double phi) const
G4double cMajorantSurfaceParameterTable[980][2]
G4double aMajorantSurfaceParameterTable[980][2]