74 const G4int arrayDim = 980;
85 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 = std::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",
116 for(
G4int i=0 ; i<arrayDim ;++i){
117 infile >>
beta >> aRead >> cRead;
144 G4double beta = std::sqrt((gamma - 1)*(gamma + 1))/gamma;
184 G4double crossSectionMajorantFunctionValue = 0;
197 theta=std::sqrt(((
G4Exp(rand2*std::log(1+cBeta*
pi*
pi)))-1)/cBeta);
198 crossSectionMajorantFunctionValue =
203 theta = std::sqrt(((
G4Exp(rand2*std::log(1+cBeta*
pi*
pi)))-1)/cBeta);
204 crossSectionMajorantFunctionValue =
209 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
210 if(crossSectionValue < 0.0) { crossSectionValue = maxBeta; }
212 }
while(maxBeta > crossSectionValue || theta >
CLHEP::pi);
225 G4double crossSectionMajorantFunctionValue = 0;
226 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
227 return crossSectionMajorantFunctionValue;
240 G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
241 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
242 G4double cosTheta = std::cos(theta);
243 G4double sinTheta2 = std::sin(theta)*std::sin(theta);
244 G4double cosPhi2 = std::cos(phi)*std::cos(phi);
250 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
251 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
252 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
254 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
255 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*
beta/oneBeta2 * cosTheta * cosPhi2
256 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
257 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*
beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
258 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (
beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
259 (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);
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);
314 if(!(polarization.
isOrthogonal(direction,kTolerance)) || mS == 0){
325 polarization2 = c.
unit();
326 mS = polarization2.
mag();
331 polarization2 = polarization
332 - polarization.
dot(direction)/direction.
dot(direction) * direction;
337 polarization2 = polarization2/mS;
356 if(shellId > 0) { level = 1; }
381 else if(k == indexMax)
388 *majorantSurfaceParameterA = aBeta;
389 *majorantSurfaceParameterC = cBeta;
404 G4ThreeVector outgoingDirection = rotation*samplingDirection;
405 return outgoingDirection;
413 G4cout <<
"Polarized Photoelectric Angular Generator" <<
G4endl;
414 G4cout <<
"PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" <<
G4endl;
415 G4cout <<
"Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." <<
G4endl;
416 G4cout <<
"For higher shells the L1 cross-section is used." <<
G4endl;
417 G4cout <<
"(see Physics Reference Manual) \n" <<
G4endl;
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
static constexpr double twopi
static constexpr double pi
static const G4double angle[DIMMOTT]
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
Hep3Vector cross(const Hep3Vector &) const
double dot(const Hep3Vector &) const
bool isOrthogonal(const Hep3Vector &v, double epsilon=tolerance) const
double howOrthogonal(const Hep3Vector &v) const
const G4ThreeVector & GetMomentumDirection() const
const G4ThreeVector & GetPolarization() const
G4double DSigmaL1shellGavrila(G4double beta, G4double theta, G4double phi) const
void PhotoElectronGetMajorantSurfaceAandCParameters(G4int shellId, G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
G4ThreeVector PerpendicularVector(const G4ThreeVector &a) const
G4ThreeVector PhotoElectronComputeFinalDirection(const G4RotationMatrix &rotation, G4double theta, G4double phi) const
void PrintGeneratorInformation() const override
G4double CrossSectionMajorantFunction(G4double theta, G4double cBeta) const
void PhotoElectronGeneratePhiAndTheta(G4int shellId, G4double beta, G4double aBeta, G4double cBeta, G4double *pphi, G4double *ptheta) const
G4double cMajorantSurfaceParameterTable[980][2]
G4RotationMatrix PhotoElectronRotationMatrix(const G4ThreeVector &direction, const G4ThreeVector &polarization)
G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double eKinEnergy, G4int shellId, const G4Material *mat=nullptr) override
~G4PhotoElectricAngularGeneratorPolarized()
G4double aMajorantSurfaceParameterTable[980][2]
G4PhotoElectricAngularGeneratorPolarized()
G4double DSigmaKshellGavrila1959(G4double beta, G4double theta, G4double phi) const
G4ThreeVector fLocalDirection
static constexpr double pi
T max(const T t1, const T t2)
brief Return the largest of the two arguments