Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // 26 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4PhotoElectricAngularGenera 32 // File name: G4PhotoElectricAngularGeneratorPolarized 33 // 33 // 34 // Author: A. C. Farinha, L. Peralta, P. Rodri 34 // Author: A. C. Farinha, L. Peralta, P. Rodrigues and A. Trindade 35 // 35 // 36 // Creation date: 36 // Creation date: 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 10 January 2006 39 // 10 January 2006 40 // 06 May 2011, Replace FILE with std::ifstrea 40 // 06 May 2011, Replace FILE with std::ifstream 41 // 41 // 42 // Class Description: 42 // Class Description: 43 // 43 // 44 // Concrete class for PhotoElectric Electron A 44 // Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation 45 // 45 // 46 // Class Description: 46 // Class Description: 47 // PhotoElectric Electron Angular Generator ba 47 // PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution. 48 // Includes polarization effects for K and L1 48 // Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961). 49 // For higher shells the L1 cross-section is u 49 // For higher shells the L1 cross-section is used. 50 // 50 // 51 // The Gavrila photoelectron angular distribut 51 // The Gavrila photoelectron angular distribution is a complex function which can not be sampled using 52 // the inverse-transform method (James 1980). 52 // the inverse-transform method (James 1980). Instead a more general approach based on the one already 53 // used to sample bremsstrahlung 2BN cross sec 53 // used to sample bremsstrahlung 2BN cross section (G4Generator2BN, Peralta, 2005) was used. 54 // 54 // 55 // M. Gavrila, "Relativistic K-Shell Photoeffe 55 // M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526 (1959) 56 // M. Gavrila, "Relativistic L-Shell Photoeffe 56 // M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961) 57 // F. James, Rept. on Prog. in Phys. 43, 1145 57 // F. James, Rept. on Prog. in Phys. 43, 1145 (1980) 58 // L. Peralta et al., "A new low-energy bremss 58 // L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005) 59 // 59 // 60 // 60 // 61 // ------------------------------------------- 61 // ------------------------------------------------------------------- 62 // 62 // 63 // 63 // 64 64 65 #include "G4PhotoElectricAngularGeneratorPolar 65 #include "G4PhotoElectricAngularGeneratorPolarized.hh" 66 #include "G4PhysicalConstants.hh" 66 #include "G4PhysicalConstants.hh" 67 #include "G4RotationMatrix.hh" 67 #include "G4RotationMatrix.hh" 68 #include "Randomize.hh" 68 #include "Randomize.hh" 69 #include "G4Exp.hh" 69 #include "G4Exp.hh" 70 70 71 G4PhotoElectricAngularGeneratorPolarized::G4Ph 71 G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized() 72 :G4VEmAngularDistribution("AngularGenSauterG 72 :G4VEmAngularDistribution("AngularGenSauterGavrilaPolarized") 73 { 73 { 74 const G4int arrayDim = 980; 74 const G4int arrayDim = 980; 75 75 76 //minimum electron beta parameter allowed 76 //minimum electron beta parameter allowed 77 betaArray[0] = 0.02; 77 betaArray[0] = 0.02; 78 //beta step 78 //beta step 79 betaArray[1] = 0.001; 79 betaArray[1] = 0.001; 80 //maximum index array for a and c tables 80 //maximum index array for a and c tables 81 betaArray[2] = arrayDim - 1; 81 betaArray[2] = arrayDim - 1; 82 82 83 // read Majorant Surface Parameters. This ar 83 // read Majorant Surface Parameters. This are required in order to generate 84 //Gavrila angular photoelectron distribution 84 //Gavrila angular photoelectron distribution 85 for(G4int level = 0; level < 2; level++){ 85 for(G4int level = 0; level < 2; level++){ 86 char nameChar0[100] = "ftab0.dat"; // K-sh 86 char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters 87 char nameChar1[100] = "ftab1.dat"; // L-sh 87 char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters 88 88 89 G4String filename; 89 G4String filename; 90 if(level == 0) filename = nameChar0; 90 if(level == 0) filename = nameChar0; 91 if(level == 1) filename = nameChar1; 91 if(level == 1) filename = nameChar1; 92 92 93 const char* path = G4FindDataDir("G4LEDATA << 93 char* path = std::getenv("G4LEDATA"); 94 if (!path) 94 if (!path) 95 { 95 { 96 G4String excep = "G4EMDataSet - G4LEDA 96 G4String excep = "G4EMDataSet - G4LEDATA environment variable not set"; 97 G4Exception("G4PhotoElectricAngularGen 97 G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized", 98 "em0006",FatalException,"G4LEDATA envi 98 "em0006",FatalException,"G4LEDATA environment variable not set"); 99 return; 99 return; 100 } 100 } 101 101 102 G4String pathString(path); 102 G4String pathString(path); 103 G4String dirFile = pathString + "/photoele 103 G4String dirFile = pathString + "/photoelectric_angular/" + filename; 104 std::ifstream infile(dirFile); 104 std::ifstream infile(dirFile); 105 if (!infile.is_open()) 105 if (!infile.is_open()) 106 { 106 { 107 G4String excep = "data file: " + dirFile + " 107 G4String excep = "data file: " + dirFile + " not found"; 108 G4Exception("G4PhotoElectricAngularGen 108 G4Exception("G4PhotoElectricAngularGeneratorPolarized::G4PhotoElectricAngularGeneratorPolarized", 109 "em0003",FatalException,excep); 109 "em0003",FatalException,excep); 110 return; 110 return; 111 } 111 } 112 112 113 // Read parameters into tables. The parame 113 // Read parameters into tables. The parameters are function of incident electron 114 // energy and shell level 114 // energy and shell level 115 G4float aRead=0,cRead=0, beta=0; 115 G4float aRead=0,cRead=0, beta=0; 116 for(G4int i=0 ; i<arrayDim ;++i){ 116 for(G4int i=0 ; i<arrayDim ;++i){ 117 infile >> beta >> aRead >> cRead; 117 infile >> beta >> aRead >> cRead; 118 aMajorantSurfaceParameterTable[i][level] 118 aMajorantSurfaceParameterTable[i][level] = aRead; 119 cMajorantSurfaceParameterTable[i][level] 119 cMajorantSurfaceParameterTable[i][level] = cRead; 120 } 120 } 121 infile.close(); 121 infile.close(); 122 } 122 } 123 } 123 } 124 124 125 //....oooOO0OOooo........oooOO0OOooo........oo 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 126 126 127 G4PhotoElectricAngularGeneratorPolarized::~G4P 127 G4PhotoElectricAngularGeneratorPolarized::~G4PhotoElectricAngularGeneratorPolarized() 128 {} 128 {} 129 129 130 //....oooOO0OOooo........oooOO0OOooo........oo 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 131 131 132 G4ThreeVector& 132 G4ThreeVector& 133 G4PhotoElectricAngularGeneratorPolarized::Samp 133 G4PhotoElectricAngularGeneratorPolarized::SampleDirection( 134 const G4Dynam 134 const G4DynamicParticle* dp, 135 G4double eKinEnergy, 135 G4double eKinEnergy, 136 G4int shellId, 136 G4int shellId, 137 const G4Material*) 137 const G4Material*) 138 { 138 { 139 // (shellId == 0) - K-shell - Polarized mod 139 // (shellId == 0) - K-shell - Polarized model for K-shell 140 // (shellId > 0) - L1-shell - Polarized mod 140 // (shellId > 0) - L1-shell - Polarized model for L1 and higher shells 141 141 142 // Calculate Lorentz term (gamma) and beta p 142 // Calculate Lorentz term (gamma) and beta parameters 143 G4double gamma = 1 + eKinEnergy/electron_mas 143 G4double gamma = 1 + eKinEnergy/electron_mass_c2; 144 G4double beta = std::sqrt((gamma - 1)*(gamm 144 G4double beta = std::sqrt((gamma - 1)*(gamma + 1))/gamma; 145 145 146 const G4ThreeVector& direction = dp->GetMome 146 const G4ThreeVector& direction = dp->GetMomentumDirection(); 147 const G4ThreeVector& polarization = dp->GetP 147 const G4ThreeVector& polarization = dp->GetPolarization(); 148 148 149 G4double theta, phi = 0; 149 G4double theta, phi = 0; 150 // Majorant surface parameters 150 // Majorant surface parameters 151 // function of the outgoing electron kinetic 151 // function of the outgoing electron kinetic energy 152 G4double aBeta = 0; 152 G4double aBeta = 0; 153 G4double cBeta = 0; 153 G4double cBeta = 0; 154 154 155 // For the outgoing kinetic energy 155 // For the outgoing kinetic energy 156 // find the current majorant surface paramet 156 // find the current majorant surface parameters 157 PhotoElectronGetMajorantSurfaceAandCParamete 157 PhotoElectronGetMajorantSurfaceAandCParameters(shellId, beta, &aBeta, &cBeta); 158 158 159 // Generate pho and theta according to the s 159 // Generate pho and theta according to the shell level 160 // and beta parameter of the electron 160 // and beta parameter of the electron 161 PhotoElectronGeneratePhiAndTheta(shellId, be 161 PhotoElectronGeneratePhiAndTheta(shellId, beta, aBeta, cBeta, &phi, &theta); 162 162 163 // Determine the rotation matrix 163 // Determine the rotation matrix 164 const G4RotationMatrix rotation = 164 const G4RotationMatrix rotation = 165 PhotoElectronRotationMatrix(direction, pol 165 PhotoElectronRotationMatrix(direction, polarization); 166 166 167 // Compute final direction of the outgoing e 167 // Compute final direction of the outgoing electron 168 fLocalDirection = PhotoElectronComputeFinalD 168 fLocalDirection = PhotoElectronComputeFinalDirection(rotation, theta, phi); 169 169 170 return fLocalDirection; 170 return fLocalDirection; 171 } 171 } 172 172 173 //....oooOO0OOooo........oooOO0OOooo........oo 173 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 174 174 175 void 175 void 176 G4PhotoElectricAngularGeneratorPolarized::Phot 176 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta( 177 G4int shellLevel, G4double beta, G4doubl 177 G4int shellLevel, G4double beta, G4double aBeta, G4double cBeta, 178 G4double *pphi, G4double *ptheta) const 178 G4double *pphi, G4double *ptheta) const 179 { 179 { 180 G4double rand1, rand2, rand3 = 0; 180 G4double rand1, rand2, rand3 = 0; 181 G4double phi = 0; 181 G4double phi = 0; 182 G4double theta = 0; 182 G4double theta = 0; 183 G4double crossSectionValue = 0; 183 G4double crossSectionValue = 0; 184 G4double crossSectionMajorantFunctionValue = 184 G4double crossSectionMajorantFunctionValue = 0; 185 G4double maxBeta = 0; 185 G4double maxBeta = 0; 186 186 187 do { 187 do { 188 188 189 rand1 = G4UniformRand(); 189 rand1 = G4UniformRand(); 190 rand2 = G4UniformRand(); 190 rand2 = G4UniformRand(); 191 rand3 = G4UniformRand(); 191 rand3 = G4UniformRand(); 192 192 193 phi=2*pi*rand1; 193 phi=2*pi*rand1; 194 194 195 if(shellLevel == 0){ 195 if(shellLevel == 0){ 196 // Polarized Gavrila Cross-Section for K 196 // Polarized Gavrila Cross-Section for K-shell (1959) 197 theta=std::sqrt(((G4Exp(rand2*std::log(1 197 theta=std::sqrt(((G4Exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta); 198 crossSectionMajorantFunctionValue = 198 crossSectionMajorantFunctionValue = 199 CrossSectionMajorantFunction(theta, cBeta); 199 CrossSectionMajorantFunction(theta, cBeta); 200 crossSectionValue = DSigmaKshellGavrila1 200 crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi); 201 } else { 201 } else { 202 // Polarized Gavrila Cross-Section for 202 // Polarized Gavrila Cross-Section for other shells (L1-shell) (1961) 203 theta = std::sqrt(((G4Exp(rand2*std::log 203 theta = std::sqrt(((G4Exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta); 204 crossSectionMajorantFunctionValue = 204 crossSectionMajorantFunctionValue = 205 CrossSectionMajorantFunction(theta, cBeta); 205 CrossSectionMajorantFunction(theta, cBeta); 206 crossSectionValue = DSigmaL1shellGavrila 206 crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi); 207 } 207 } 208 208 209 maxBeta=rand3*aBeta*crossSectionMajorantFu 209 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue; 210 if(crossSectionValue < 0.0) { crossSection 210 if(crossSectionValue < 0.0) { crossSectionValue = maxBeta; } 211 211 212 } while(maxBeta > crossSectionValue || theta 212 } while(maxBeta > crossSectionValue || theta > CLHEP::pi); 213 213 214 *pphi = phi; 214 *pphi = phi; 215 *ptheta = theta; 215 *ptheta = theta; 216 } 216 } 217 217 218 //....oooOO0OOooo........oooOO0OOooo........oo 218 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 219 219 220 G4double 220 G4double 221 G4PhotoElectricAngularGeneratorPolarized::Cros 221 G4PhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction( 222 G4double theta, G4double cBeta) const 222 G4double theta, G4double cBeta) const 223 { 223 { 224 // Compute Majorant Function 224 // Compute Majorant Function 225 G4double crossSectionMajorantFunctionValue = 225 G4double crossSectionMajorantFunctionValue = 0; 226 crossSectionMajorantFunctionValue = theta/(1 226 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta); 227 return crossSectionMajorantFunctionValue; 227 return crossSectionMajorantFunctionValue; 228 } 228 } 229 229 230 //....oooOO0OOooo........oooOO0OOooo........oo 230 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 231 231 232 G4double 232 G4double 233 G4PhotoElectricAngularGeneratorPolarized::DSig 233 G4PhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959( 234 G4double beta, G4double theta, G4doub 234 G4double beta, G4double theta, G4double phi) const 235 { 235 { 236 //Double differential K shell cross-section 236 //Double differential K shell cross-section (Gavrila 1959) 237 237 238 G4double beta2 = beta*beta; 238 G4double beta2 = beta*beta; 239 G4double oneBeta2 = 1 - beta2; 239 G4double oneBeta2 = 1 - beta2; 240 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 240 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 241 G4double oneBeta2_to_3_2 = std::pow(oneBeta2 241 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5); 242 G4double cosTheta = std::cos(theta); 242 G4double cosTheta = std::cos(theta); 243 G4double sinTheta2 = std::sin(theta)*std::si 243 G4double sinTheta2 = std::sin(theta)*std::sin(theta); 244 G4double cosPhi2 = std::cos(phi)*std::cos(ph 244 G4double cosPhi2 = std::cos(phi)*std::cos(phi); 245 G4double oneBetaCosTheta = 1-beta*cosTheta; 245 G4double oneBetaCosTheta = 1-beta*cosTheta; 246 G4double dsigma = 0; 246 G4double dsigma = 0; 247 G4double firstTerm = 0; 247 G4double firstTerm = 0; 248 G4double secondTerm = 0; 248 G4double secondTerm = 0; 249 249 250 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe 250 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) * 251 (sinTheta2 * cosPhi2)/std::pow(o 251 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)* 252 (1-sqrtOneBeta2)/(4*oneBeta2_to_ 252 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3); 253 253 254 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st 254 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) * 255 (4*beta2/sqrtOneBeta2 * sinThet 255 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2 256 - 4*(1-sqrtOneBeta2)/oneBeta2 * 256 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta 257 + 4*beta2*(1-sqrtOneBeta2)/oneB 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*one 258 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 + 259 (1-sqrtOneBeta2)/oneBeta2_to_3_ 259 (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2); 260 260 261 dsigma = ( firstTerm*(1-pi*fine_structure_co 261 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) )*std::sin(theta); 262 262 263 return dsigma; 263 return dsigma; 264 } 264 } 265 265 266 //....oooOO0OOooo........oooOO0OOooo........oo 266 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 267 267 268 G4double 268 G4double 269 G4PhotoElectricAngularGeneratorPolarized::DSig 269 G4PhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila( 270 G4double beta, G4double theta, G4doubl 270 G4double beta, G4double theta, G4double phi) const 271 { 271 { 272 //Double differential L1 shell cross-section 272 //Double differential L1 shell cross-section (Gavrila 1961) 273 // 26Oct2022: included factor (1/8) in dsigm << 274 G4double beta2 = beta*beta; 273 G4double beta2 = beta*beta; 275 G4double oneBeta2 = 1-beta2; 274 G4double oneBeta2 = 1-beta2; 276 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 275 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 277 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1 276 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5); 278 G4double cosTheta = std::cos(theta); 277 G4double cosTheta = std::cos(theta); 279 G4double sinTheta2 =std::sin(theta)*std::sin 278 G4double sinTheta2 =std::sin(theta)*std::sin(theta); 280 G4double cosPhi2 = std::cos(phi)*std::cos(ph 279 G4double cosPhi2 = std::cos(phi)*std::cos(phi); 281 G4double oneBetaCosTheta = 1-beta*cosTheta; 280 G4double oneBetaCosTheta = 1-beta*cosTheta; 282 281 283 G4double dsigma = 0; 282 G4double dsigma = 0; 284 G4double firstTerm = 0; 283 G4double firstTerm = 0; 285 G4double secondTerm = 0; 284 G4double secondTerm = 0; 286 285 287 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe 286 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) 288 * (sinTheta2 * cosPhi2)/std::po 287 * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)* 289 (1-sqrtOneBeta2)/(4*oneBeta2_to_ 288 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3); 290 289 291 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st 290 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) * 292 (4*beta2/sqrtOneBeta2 * sinThet 291 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2 293 - 4*(1-sqrtOneBeta2)/oneBeta2 * 292 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta 294 + 4*beta2*(1-sqrtOneBeta2)/oneB 293 + 4*beta2*(1-sqrtOneBeta2)/oneBeta2_to_3_2 - 4*beta*(1-sqrtOneBeta2)*(1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta) 295 + (1-sqrtOneBeta2)/(4*beta2*one 294 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 + 296 (1-sqrtOneBeta2)/oneBeta2_to_3_ 295 (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2); 297 296 298 dsigma = ( firstTerm*(1-pi*fine_structure_co << 297 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) )*std::sin(theta); 299 298 300 return dsigma; 299 return dsigma; 301 } 300 } 302 301 303 //....oooOO0OOooo........oooOO0OOooo........oo 302 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 304 303 305 G4RotationMatrix 304 G4RotationMatrix 306 G4PhotoElectricAngularGeneratorPolarized::Phot 305 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix( 307 const G4ThreeVector& direction, 306 const G4ThreeVector& direction, 308 const G4ThreeVector& polarization) 307 const G4ThreeVector& polarization) 309 { 308 { 310 G4double mK = direction.mag(); 309 G4double mK = direction.mag(); 311 G4double mS = polarization.mag(); 310 G4double mS = polarization.mag(); 312 G4ThreeVector polarization2 = polarization; 311 G4ThreeVector polarization2 = polarization; 313 const G4double kTolerance = 1e-6; 312 const G4double kTolerance = 1e-6; 314 313 315 if(!(polarization.isOrthogonal(direction,kTo 314 if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){ 316 G4ThreeVector d0 = direction.unit(); 315 G4ThreeVector d0 = direction.unit(); 317 G4ThreeVector a1 = PerpendicularVector(d0) 316 G4ThreeVector a1 = PerpendicularVector(d0); 318 G4ThreeVector a0 = a1.unit(); 317 G4ThreeVector a0 = a1.unit(); 319 G4double rand1 = G4UniformRand(); 318 G4double rand1 = G4UniformRand(); 320 G4double angle = twopi*rand1; 319 G4double angle = twopi*rand1; 321 G4ThreeVector b0 = d0.cross(a0); 320 G4ThreeVector b0 = d0.cross(a0); 322 G4ThreeVector c; 321 G4ThreeVector c; 323 c.setX(std::cos(angle)*(a0.x())+std::sin(a 322 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x()); 324 c.setY(std::cos(angle)*(a0.y())+std::sin(a 323 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y()); 325 c.setZ(std::cos(angle)*(a0.z())+std::sin(a 324 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z()); 326 polarization2 = c.unit(); 325 polarization2 = c.unit(); 327 mS = polarization2.mag(); 326 mS = polarization2.mag(); 328 }else 327 }else 329 { 328 { 330 if ( polarization.howOrthogonal(directio 329 if ( polarization.howOrthogonal(direction) != 0) 331 { 330 { 332 polarization2 = polarization 331 polarization2 = polarization 333 - polarization.dot(direction)/direction. 332 - polarization.dot(direction)/direction.dot(direction) * direction; 334 } 333 } 335 } 334 } 336 335 337 G4ThreeVector direction2 = direction/mK; 336 G4ThreeVector direction2 = direction/mK; 338 polarization2 = polarization2/mS; 337 polarization2 = polarization2/mS; 339 338 340 G4ThreeVector y = direction2.cross(polarizat 339 G4ThreeVector y = direction2.cross(polarization2); 341 340 342 G4RotationMatrix R(polarization2,y,direction 341 G4RotationMatrix R(polarization2,y,direction2); 343 return R; 342 return R; 344 } 343 } 345 344 346 //....oooOO0OOooo........oooOO0OOooo........oo 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 347 346 348 void 347 void 349 G4PhotoElectricAngularGeneratorPolarized::Phot 348 G4PhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(G4int shellId, G4double beta, G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const 350 { 349 { 351 // This member function finds for a given sh 350 // This member function finds for a given shell and beta value 352 // of the outgoing electron the correct Majo 351 // of the outgoing electron the correct Majorant Surface parameters 353 G4double aBeta,cBeta; 352 G4double aBeta,cBeta; 354 G4double bMin,bStep; 353 G4double bMin,bStep; 355 G4int indexMax; 354 G4int indexMax; 356 G4int level = 0; 355 G4int level = 0; 357 if(shellId > 0) { level = 1; } 356 if(shellId > 0) { level = 1; } 358 357 359 bMin = betaArray[0]; 358 bMin = betaArray[0]; 360 bStep = betaArray[1]; 359 bStep = betaArray[1]; 361 indexMax = (G4int)betaArray[2]; 360 indexMax = (G4int)betaArray[2]; 362 const G4double kBias = 1e-9; 361 const G4double kBias = 1e-9; 363 362 364 G4int k = (G4int)((beta-bMin+kBias)/bStep); 363 G4int k = (G4int)((beta-bMin+kBias)/bStep); 365 364 366 if(k < 0) 365 if(k < 0) 367 k = 0; 366 k = 0; 368 if(k > indexMax) 367 if(k > indexMax) 369 k = indexMax; 368 k = indexMax; 370 369 371 if(k == 0) 370 if(k == 0) 372 aBeta = std::max(aMajorantSurfaceParameter 371 aBeta = std::max(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]); 373 else if(k==indexMax) 372 else if(k==indexMax) 374 aBeta = std::max(aMajorantSurfaceParameter 373 aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]); 375 else{ 374 else{ 376 aBeta = std::max(aMajorantSurfaceParameter 375 aBeta = std::max(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]); 377 aBeta = std::max(aBeta,aMajorantSurfacePar 376 aBeta = std::max(aBeta,aMajorantSurfaceParameterTable[k+1][level]); 378 } 377 } 379 378 380 if(k == 0) 379 if(k == 0) 381 cBeta = std::max(cMajorantSurfaceParameter 380 cBeta = std::max(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]); 382 else if(k == indexMax) 381 else if(k == indexMax) 383 cBeta = std::max(cMajorantSurfaceParameter 382 cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]); 384 else{ 383 else{ 385 cBeta = std::max(cMajorantSurfaceParameter 384 cBeta = std::max(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]); 386 cBeta = std::max(cBeta,cMajorantSurfacePar 385 cBeta = std::max(cBeta,cMajorantSurfaceParameterTable[k+1][level]); 387 } 386 } 388 387 389 *majorantSurfaceParameterA = aBeta; 388 *majorantSurfaceParameterA = aBeta; 390 *majorantSurfaceParameterC = cBeta; 389 *majorantSurfaceParameterC = cBeta; 391 } 390 } 392 391 393 //....oooOO0OOooo........oooOO0OOooo........oo 392 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 394 393 395 G4ThreeVector G4PhotoElectricAngularGeneratorP 394 G4ThreeVector G4PhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, G4double theta, G4double phi) const 396 { 395 { 397 //computes the photoelectron momentum unitar 396 //computes the photoelectron momentum unitary vector 398 G4double sint = std::sin(theta); 397 G4double sint = std::sin(theta); 399 G4double px = std::cos(phi)*sint; 398 G4double px = std::cos(phi)*sint; 400 G4double py = std::sin(phi)*sint; 399 G4double py = std::sin(phi)*sint; 401 G4double pz = std::cos(theta); 400 G4double pz = std::cos(theta); 402 401 403 G4ThreeVector samplingDirection(px,py,pz); 402 G4ThreeVector samplingDirection(px,py,pz); 404 403 405 G4ThreeVector outgoingDirection = rotation*s 404 G4ThreeVector outgoingDirection = rotation*samplingDirection; 406 return outgoingDirection; 405 return outgoingDirection; 407 } 406 } 408 407 409 //....oooOO0OOooo........oooOO0OOooo........oo 408 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 410 409 411 void G4PhotoElectricAngularGeneratorPolarized: 410 void G4PhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const 412 { 411 { 413 G4cout << "\n" << G4endl; 412 G4cout << "\n" << G4endl; 414 G4cout << "Polarized Photoelectric Angular G 413 G4cout << "Polarized Photoelectric Angular Generator" << G4endl; 415 G4cout << "PhotoElectric Electron Angular Ge 414 G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl; 416 G4cout << "Includes polarization effects for 415 G4cout << "Includes polarization effects for K and L1 atomic shells, according to Gavrilla (1959, 1961)." << G4endl; 417 G4cout << "For higher shells the L1 cross-se 416 G4cout << "For higher shells the L1 cross-section is used." << G4endl; 418 G4cout << "(see Physics Reference Manual) \n 417 G4cout << "(see Physics Reference Manual) \n" << G4endl; 419 } 418 } 420 419 421 //....oooOO0OOooo........oooOO0OOooo........oo 420 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 422 421 423 G4ThreeVector 422 G4ThreeVector 424 G4PhotoElectricAngularGeneratorPolarized::Perp 423 G4PhotoElectricAngularGeneratorPolarized::PerpendicularVector( 425 const G4ThreeVector& a) const 424 const G4ThreeVector& a) const 426 { 425 { 427 G4double dx = a.x(); 426 G4double dx = a.x(); 428 G4double dy = a.y(); 427 G4double dy = a.y(); 429 G4double dz = a.z(); 428 G4double dz = a.z(); 430 G4double x = dx < 0.0 ? -dx : dx; 429 G4double x = dx < 0.0 ? -dx : dx; 431 G4double y = dy < 0.0 ? -dy : dy; 430 G4double y = dy < 0.0 ? -dy : dy; 432 G4double z = dz < 0.0 ? -dz : dz; 431 G4double z = dz < 0.0 ? -dz : dz; 433 if (x < y) { 432 if (x < y) { 434 return x < z ? G4ThreeVector(-dy,dx,0) : G 433 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); 435 }else{ 434 }else{ 436 return y < z ? G4ThreeVector(dz,0,-dx) : G 435 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); 437 } 436 } 438 } 437 } 439 438