physics/src/G4RDPhotoElectricAngularGeneratorPolarized.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // -------------------------------------------------------------------
0028 //
0029 // GEANT4 Class file
0030 //
0031 //
0032 // File name:     G4RDPhotoElectricAngularGeneratorPolarized
0033 //
0034 // Author: A. C. Farinha, L. Peralta, P. Rodrigues and A. Trindade
0035 //
0036 // Creation date:
0037 //
0038 // Modifications:
0039 // 10 January 2006
0040 //
0041 // Class Description:
0042 //
0043 // Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation
0044 //
0045 // Class Description:
0046 // PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution.
0047 // Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961).
0048 // For higher shells the L1 cross-section is used.
0049 //
0050 // The Gavrila photoelectron angular distribution is a complex function which can not be sampled using
0051 // the inverse-transform method (James 1980). Instead a more general approach based on the one already
0052 // used to sample bremsstrahlung 2BN cross section (G4RDGenerator2BN, Peralta, 2005) was used.
0053 //
0054 // M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526   (1959)
0055 // M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961)
0056 // F. James, Rept. on Prog. in Phys. 43, 1145 (1980)
0057 // L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005)
0058 //
0059 //
0060 // -------------------------------------------------------------------
0061 //
0062 //
0063
0064 #include "G4RDPhotoElectricAngularGeneratorPolarized.hh"
0065 #include "G4RotationMatrix.hh"
0066 #include "Randomize.hh"
0067 #include "G4PhysicalConstants.hh"
0068
0069 //
0070
0071 G4RDPhotoElectricAngularGeneratorPolarized::G4RDPhotoElectricAngularGeneratorPolarized(const G4String& name):G4RDVPhotoElectricAngularDistribution(name)
0072 {
0073   const G4int arrayDim = 980;
0074
0075   //minimum electron beta parameter allowed
0076   betaArray[0] = 0.02;
0077   //beta step
0078   betaArray[1] = 0.001;
0079   //maximum index array for a and c tables
0080   betaArray[2] = arrayDim - 1;
0081
0082   // read Majorant Surface Parameters. This are required in order to generate Gavrila angular photoelectron distribution
0083   for(G4int level = 0; level < 2; level++){
0084
0085     char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters
0086     char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters
0087
0088     G4String filename;
0089     if(level == 0) filename = nameChar0;
0090     if(level == 1) filename = nameChar1;
0091
0092     const char* path = G4FindDataDir("G4LEDATA");
0093     if (!path)
0094       {
0095         G4String excep = "G4LEDATA environment variable not set!";
0096         G4Exception("G4RDPhotoElectricAngularGeneratorPolarized()",
0097                     "InvalidSetup", FatalException, excep);
0098       }
0099
0100     G4String pathString(path);
0101     G4String dirFile = pathString + "/photoelectric_angular/" + filename;
0102     FILE *infile;
0103     infile = fopen(dirFile,"r");
0104     if (infile == 0)
0105       {
0106     G4String excep = "Data file: " + dirFile + " not found";
0107     G4Exception("G4RDPhotoElectricAngularGeneratorPolarized()",
0108                     "DataNotFound", FatalException, excep);
0109       }
0110
0111     // Read parameters into tables. The parameters are function of incident electron energy and shell level
0112     G4float aRead,cRead, beta;
0113     for(G4int i=0 ; i<arrayDim ;i++){
0114       fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead);
0115       aMajorantSurfaceParameterTable[i][level] = aRead;
0116       cMajorantSurfaceParameterTable[i][level] = cRead;
0117     }
0118     fclose(infile);
0119
0120   }
0121 }
0122
0123 //
0124
0125 G4RDPhotoElectricAngularGeneratorPolarized::~G4RDPhotoElectricAngularGeneratorPolarized()
0126 {;}
0127
0128 //
0129
0130 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::GetPhotoElectronDirection(const G4ThreeVector& direction, const G4double eKineticEnergy,
0131                                           const G4ThreeVector& polarization, const G4int shellId) const
0132 {
0133   // Calculate Lorentz term (gamma) and beta parameters
0134   G4double gamma   = 1. + eKineticEnergy/electron_mass_c2;
0135   G4double beta  = std::sqrt(gamma*gamma-1.)/gamma;
0136
0137   G4double theta, phi = 0;
0138   G4double aBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy)
0139   G4double cBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy)
0140
0141   G4int shellLevel = 0;
0142   if(shellId <  2) shellLevel = 0; // K-shell  // Polarized model for K-shell
0143   if(shellId >= 2) shellLevel = 1; // L1-shell // Polarized model for L1 and higher shells
0144
0145   // For the outgoing kinetic energy find the current majorant surface parameters
0146   PhotoElectronGetMajorantSurfaceAandCParameters( shellLevel, beta, &aBeta, &cBeta);
0147
0148   // Generate pho and theta according to the shell level and beta parameter of the electron
0149   PhotoElectronGeneratePhiAndTheta(shellLevel, beta, aBeta, cBeta, &phi, &theta);
0150
0151   // Determine the rotation matrix
0152   G4RotationMatrix rotation = PhotoElectronRotationMatrix(direction, polarization);
0153
0154   // Compute final direction of the outgoing electron
0155   G4ThreeVector final_direction = PhotoElectronComputeFinalDirection(rotation, theta, phi);
0156
0157   return final_direction;
0158 }
0159
0160 //
0161
0162 void G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(const G4int shellLevel, const G4double beta,
0163                                             const G4double aBeta, const G4double cBeta,
0164                                             G4double *pphi, G4double *ptheta) const
0165 {
0166   G4double rand1, rand2, rand3 = 0;
0167   G4double phi = 0;
0168   G4double theta = 0;
0169   G4double crossSectionValue = 0;
0170   G4double crossSectionMajorantFunctionValue = 0;
0171   G4double maxBeta = 0;
0172
0173   do {
0174
0175     rand1 = G4UniformRand();
0176     rand2 = G4UniformRand();
0177     rand3 = G4UniformRand();
0178
0179     phi=2*pi*rand1;
0180
0181     if(shellLevel == 0){
0182
0183       // Polarized Gavrila Cross-Section for K-shell (1959)
0184       theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
0185       crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta);
0186       crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi);
0187
0188     } else {
0189
0190       //  Polarized Gavrila Cross-Section for other shells (L1-shell) (1961)
0191       theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta);
0192       crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta);
0193       crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi);
0194
0195     }
0196
0197     maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue;
0198
0199   }while(maxBeta > crossSectionValue);
0200
0201   *pphi = phi;
0202   *ptheta = theta;
0203 }
0204
0205 //
0206
0207 G4double G4RDPhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(const G4double theta, const G4double cBeta) const
0208 {
0209   // Compute Majorant Function
0210   G4double crossSectionMajorantFunctionValue = 0;
0211   crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta);
0212   return crossSectionMajorantFunctionValue;
0213 }
0214
0215 //
0216
0217 G4double G4RDPhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(const G4double beta, const G4double theta, const G4double phi) const
0218 {
0219
0220   //Double differential K shell cross-section (Gavrila 1959)
0221
0222   G4double beta2 = beta*beta;
0223   G4double oneBeta2 = 1 - beta2;
0224   G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
0225   G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5);
0226   G4double cosTheta = std::cos(theta);
0227   G4double sinTheta2 = std::sin(theta)*std::sin(theta);
0228   G4double cosPhi2 = std::cos(phi)*std::cos(phi);
0229   G4double oneBetaCosTheta = 1-beta*cosTheta;
0230   G4double dsigma = 0;
0231   G4double firstTerm = 0;
0232   G4double secondTerm = 0;
0233
0234   firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) *
0235               (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
0236               (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
0237
0238   secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
0239                (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
0240                - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
0241                + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
0242                + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
0243                (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2);
0244
0245   dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
0246
0247   return dsigma;
0248 }
0249
0250 //
0251
0252 G4double G4RDPhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(const G4double beta, const G4double theta, const G4double phi) const
0253 {
0254
0255   //Double differential L1 shell cross-section (Gavrila 1961)
0256
0257   G4double beta2 = beta*beta;
0258   G4double oneBeta2 = 1-beta2;
0259   G4double sqrtOneBeta2 = std::sqrt(oneBeta2);
0260   G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5);
0261   G4double cosTheta = std::cos(theta);
0262   G4double sinTheta2 =std::sin(theta)*std::sin(theta);
0263   G4double cosPhi2 = std::cos(phi)*std::cos(phi);
0264   G4double oneBetaCosTheta = 1-beta*cosTheta;
0265
0266   G4double dsigma = 0;
0267   G4double firstTerm = 0;
0268   G4double secondTerm = 0;
0269
0270   firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2)
0271               *  (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)*
0272               (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3);
0273
0274   secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) *
0275                (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2
0276                - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta
0277                + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta)
0278                + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 +
0279                (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2);
0280
0281   dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) );
0282
0283   return dsigma;
0284 }
0285
0286 G4double G4RDPhotoElectricAngularGeneratorPolarized::GetMax(const G4double arg1, const G4double arg2) const
0287 {
0288   if (arg1 > arg2)
0289     return arg1;
0290   else
0291     return arg2;
0292 }
0293
0294 //
0295
0296 G4RotationMatrix G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(const G4ThreeVector& direction,
0297                                                const G4ThreeVector& polarization) const
0298 {
0299   G4double mK = direction.mag();
0300   G4double mS = polarization.mag();
0301   G4ThreeVector polarization2 = polarization;
0302   const G4double kTolerance = 1e-6;
0303
0304   if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){
0305     G4ThreeVector d0 = direction.unit();
0306     G4ThreeVector a1 = SetPerpendicularVector(d0);
0307     G4ThreeVector a0 = a1.unit();
0308     G4double rand1 = G4UniformRand();
0309     G4double angle = twopi*rand1;
0310     G4ThreeVector b0 = d0.cross(a0);
0311     G4ThreeVector c;
0312     c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
0313     c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
0314     c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
0315     polarization2 = c.unit();
0316     mS = polarization2.mag();
0317   }else
0318     {
0319       if ( polarization.howOrthogonal(direction) != 0)
0320     {
0321       polarization2 = polarization - polarization.dot(direction)/direction.dot(direction) * direction;
0322     }
0323     }
0324
0325   G4ThreeVector direction2 = direction/mK;
0326   polarization2 = polarization2/mS;
0327
0328   G4ThreeVector y = direction2.cross(polarization2);
0329
0330   G4RotationMatrix R(polarization2,y,direction2);
0331   return R;
0332 }
0333
0334 void G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(const G4int shellLevel, const G4double beta,G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const
0335 {
0336   // This member function finds for a given shell and beta value of the outgoing electron the correct Majorant Surface parameters
0337
0338   G4double aBeta,cBeta;
0339   G4double bMin,bStep;
0340   G4int indexMax;
0341   G4int level = shellLevel;
0342   if(shellLevel > 1) level = 1; // protection since only K and L1 polarized double differential cross-sections were implemented
0343
0344   bMin = betaArray[0];
0345   bStep = betaArray[1];
0346   indexMax = (G4int)betaArray[2];
0347   const G4double kBias = 1e-9;
0348
0349   G4int k = (G4int)((beta-bMin+kBias)/bStep);
0350
0351   if(k < 0)
0352     k = 0;
0353   if(k > indexMax)
0354     k = indexMax;
0355
0356   if(k == 0)
0357     aBeta = GetMax(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]);
0358   else if(k==indexMax)
0359     aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
0360   else{
0361     aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]);
0362     aBeta = GetMax(aBeta,aMajorantSurfaceParameterTable[k+1][level]);
0363   }
0364
0365   if(k == 0)
0366     cBeta = GetMax(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]);
0367   else if(k == indexMax)
0368     cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
0369   else{
0370     cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]);
0371     cBeta = GetMax(cBeta,cMajorantSurfaceParameterTable[k+1][level]);
0372   }
0373
0374   *majorantSurfaceParameterA = aBeta;
0375   *majorantSurfaceParameterC = cBeta;
0376
0377 }
0378
0379
0380 //
0381 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, const G4double theta, const G4double phi) const
0382 {
0383
0384   //computes the photoelectron momentum unitary vector
0385   G4double px = std::cos(phi)*std::sin(theta);
0386   G4double py = std::sin(phi)*std::sin(theta);
0387   G4double pz = std::cos(theta);
0388
0389   G4ThreeVector samplingDirection(px,py,pz);
0390
0391   G4ThreeVector outgoingDirection = rotation*samplingDirection;
0392   return outgoingDirection;
0393 }
0394
0395 //
0396
0397 void G4RDPhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const
0398 {
0399   G4cout << "\n" << G4endl;
0400   G4cout << "Polarized Photoelectric Angular Generator" << G4endl;
0401   G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl;
0402   G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl;
0403   G4cout << "For higher shells the L1 cross-section is used." << G4endl;
0404   G4cout << "(see Physics Reference Manual) \n" << G4endl;
0405 }
0406
0407 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::SetPerpendicularVector(const G4ThreeVector& a) const
0408 {
0409   G4double dx = a.x();
0410   G4double dy = a.y();
0411   G4double dz = a.z();
0412   G4double x = dx < 0.0 ? -dx : dx;
0413   G4double y = dy < 0.0 ? -dy : dy;
0414   G4double z = dz < 0.0 ? -dz : dz;
0415   if (x < y) {
0416     return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
0417   }else{
0418     return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
0419   }
0420 }