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: G4RDPhotoElectricAngularGene 32 // File name: G4RDPhotoElectricAngularGeneratorPolarized 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 // 40 // 41 // Class Description: 41 // Class Description: 42 // 42 // 43 // Concrete class for PhotoElectric Electron A 43 // Concrete class for PhotoElectric Electron Angular Polarized Distribution Generation 44 // 44 // 45 // Class Description: 45 // Class Description: 46 // PhotoElectric Electron Angular Generator ba 46 // PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution. 47 // Includes polarization effects for K and L1 47 // Includes polarization effects for K and L1 atomic shells, according to Gavrila (1959, 1961). 48 // For higher shells the L1 cross-section is u 48 // For higher shells the L1 cross-section is used. 49 // 49 // 50 // The Gavrila photoelectron angular distribut 50 // The Gavrila photoelectron angular distribution is a complex function which can not be sampled using 51 // the inverse-transform method (James 1980). 51 // the inverse-transform method (James 1980). Instead a more general approach based on the one already 52 // used to sample bremsstrahlung 2BN cross sec 52 // used to sample bremsstrahlung 2BN cross section (G4RDGenerator2BN, Peralta, 2005) was used. 53 // 53 // 54 // M. Gavrila, "Relativistic K-Shell Photoeffe 54 // M. Gavrila, "Relativistic K-Shell Photoeffect", Phys. Rev. 113, 514-526 (1959) 55 // M. Gavrila, "Relativistic L-Shell Photoeffe 55 // M. Gavrila, "Relativistic L-Shell Photoeffect", Phys. Rev. 124, 1132-1141 (1961) 56 // F. James, Rept. on Prog. in Phys. 43, 1145 56 // F. James, Rept. on Prog. in Phys. 43, 1145 (1980) 57 // L. Peralta et al., "A new low-energy bremss 57 // L. Peralta et al., "A new low-energy bremsstrahlung generator for GEANT4", Radiat. Prot. Dosimetry. 116, 59-64 (2005) 58 // 58 // 59 // 59 // 60 // ------------------------------------------- 60 // ------------------------------------------------------------------- 61 // 61 // 62 // 62 // 63 63 64 #include "G4RDPhotoElectricAngularGeneratorPol 64 #include "G4RDPhotoElectricAngularGeneratorPolarized.hh" 65 #include "G4RotationMatrix.hh" 65 #include "G4RotationMatrix.hh" 66 #include "Randomize.hh" 66 #include "Randomize.hh" 67 #include "G4PhysicalConstants.hh" 67 #include "G4PhysicalConstants.hh" 68 68 69 // 69 // 70 70 71 G4RDPhotoElectricAngularGeneratorPolarized::G4 71 G4RDPhotoElectricAngularGeneratorPolarized::G4RDPhotoElectricAngularGeneratorPolarized(const G4String& name):G4RDVPhotoElectricAngularDistribution(name) 72 { 72 { 73 const G4int arrayDim = 980; 73 const G4int arrayDim = 980; 74 74 75 //minimum electron beta parameter allowed 75 //minimum electron beta parameter allowed 76 betaArray[0] = 0.02; 76 betaArray[0] = 0.02; 77 //beta step 77 //beta step 78 betaArray[1] = 0.001; 78 betaArray[1] = 0.001; 79 //maximum index array for a and c tables 79 //maximum index array for a and c tables 80 betaArray[2] = arrayDim - 1; 80 betaArray[2] = arrayDim - 1; 81 81 82 // read Majorant Surface Parameters. This ar 82 // read Majorant Surface Parameters. This are required in order to generate Gavrila angular photoelectron distribution 83 for(G4int level = 0; level < 2; level++){ 83 for(G4int level = 0; level < 2; level++){ 84 84 85 char nameChar0[100] = "ftab0.dat"; // K-sh 85 char nameChar0[100] = "ftab0.dat"; // K-shell Majorant Surface Parameters 86 char nameChar1[100] = "ftab1.dat"; // L-sh 86 char nameChar1[100] = "ftab1.dat"; // L-shell Majorant Surface Parameters 87 87 88 G4String filename; 88 G4String filename; 89 if(level == 0) filename = nameChar0; 89 if(level == 0) filename = nameChar0; 90 if(level == 1) filename = nameChar1; 90 if(level == 1) filename = nameChar1; 91 91 92 const char* path = G4FindDataDir("G4LEDATA << 92 char* path = getenv("G4LEDATA"); 93 if (!path) 93 if (!path) 94 { 94 { 95 G4String excep = "G4LEDATA environment 95 G4String excep = "G4LEDATA environment variable not set!"; 96 G4Exception("G4RDPhotoElectricAngularG 96 G4Exception("G4RDPhotoElectricAngularGeneratorPolarized()", 97 "InvalidSetup", FatalExcep 97 "InvalidSetup", FatalException, excep); 98 } 98 } 99 99 100 G4String pathString(path); 100 G4String pathString(path); 101 G4String dirFile = pathString + "/photoele 101 G4String dirFile = pathString + "/photoelectric_angular/" + filename; 102 FILE *infile; 102 FILE *infile; 103 infile = fopen(dirFile,"r"); 103 infile = fopen(dirFile,"r"); 104 if (infile == 0) 104 if (infile == 0) 105 { 105 { 106 G4String excep = "Data file: " + dirFile + " 106 G4String excep = "Data file: " + dirFile + " not found"; 107 G4Exception("G4RDPhotoElectricAngularGenerat 107 G4Exception("G4RDPhotoElectricAngularGeneratorPolarized()", 108 "DataNotFound", FatalExcep 108 "DataNotFound", FatalException, excep); 109 } 109 } 110 110 111 // Read parameters into tables. The parame 111 // Read parameters into tables. The parameters are function of incident electron energy and shell level 112 G4float aRead,cRead, beta; 112 G4float aRead,cRead, beta; 113 for(G4int i=0 ; i<arrayDim ;i++){ 113 for(G4int i=0 ; i<arrayDim ;i++){ 114 fscanf(infile,"%f\t %e\t %e",&beta,&aRea 114 fscanf(infile,"%f\t %e\t %e",&beta,&aRead,&cRead); 115 aMajorantSurfaceParameterTable[i][level] 115 aMajorantSurfaceParameterTable[i][level] = aRead; 116 cMajorantSurfaceParameterTable[i][level] 116 cMajorantSurfaceParameterTable[i][level] = cRead; 117 } 117 } 118 fclose(infile); 118 fclose(infile); 119 119 120 } 120 } 121 } 121 } 122 122 123 // 123 // 124 124 125 G4RDPhotoElectricAngularGeneratorPolarized::~G 125 G4RDPhotoElectricAngularGeneratorPolarized::~G4RDPhotoElectricAngularGeneratorPolarized() 126 {;} 126 {;} 127 127 128 // 128 // 129 129 130 G4ThreeVector G4RDPhotoElectricAngularGenerato 130 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::GetPhotoElectronDirection(const G4ThreeVector& direction, const G4double eKineticEnergy, 131 const G4ThreeVector& pol 131 const G4ThreeVector& polarization, const G4int shellId) const 132 { 132 { 133 // Calculate Lorentz term (gamma) and beta p 133 // Calculate Lorentz term (gamma) and beta parameters 134 G4double gamma = 1. + eKineticEnergy/elect 134 G4double gamma = 1. + eKineticEnergy/electron_mass_c2; 135 G4double beta = std::sqrt(gamma*gamma-1.)/g 135 G4double beta = std::sqrt(gamma*gamma-1.)/gamma; 136 136 137 G4double theta, phi = 0; 137 G4double theta, phi = 0; 138 G4double aBeta = 0; // Majorant surface para 138 G4double aBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy) 139 G4double cBeta = 0; // Majorant surface para 139 G4double cBeta = 0; // Majorant surface parameter (function of the outgoing electron kinetic energy) 140 140 141 G4int shellLevel = 0; 141 G4int shellLevel = 0; 142 if(shellId < 2) shellLevel = 0; // K-shell 142 if(shellId < 2) shellLevel = 0; // K-shell // Polarized model for K-shell 143 if(shellId >= 2) shellLevel = 1; // L1-shell 143 if(shellId >= 2) shellLevel = 1; // L1-shell // Polarized model for L1 and higher shells 144 144 145 // For the outgoing kinetic energy find the 145 // For the outgoing kinetic energy find the current majorant surface parameters 146 PhotoElectronGetMajorantSurfaceAandCParamete 146 PhotoElectronGetMajorantSurfaceAandCParameters( shellLevel, beta, &aBeta, &cBeta); 147 147 148 // Generate pho and theta according to the s 148 // Generate pho and theta according to the shell level and beta parameter of the electron 149 PhotoElectronGeneratePhiAndTheta(shellLevel, 149 PhotoElectronGeneratePhiAndTheta(shellLevel, beta, aBeta, cBeta, &phi, &theta); 150 150 151 // Determine the rotation matrix 151 // Determine the rotation matrix 152 G4RotationMatrix rotation = PhotoElectronRot 152 G4RotationMatrix rotation = PhotoElectronRotationMatrix(direction, polarization); 153 153 154 // Compute final direction of the outgoing e 154 // Compute final direction of the outgoing electron 155 G4ThreeVector final_direction = PhotoElectro 155 G4ThreeVector final_direction = PhotoElectronComputeFinalDirection(rotation, theta, phi); 156 156 157 return final_direction; 157 return final_direction; 158 } 158 } 159 159 160 // 160 // 161 161 162 void G4RDPhotoElectricAngularGeneratorPolarize 162 void G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronGeneratePhiAndTheta(const G4int shellLevel, const G4double beta, 163 const G4double aBeta 163 const G4double aBeta, const G4double cBeta, 164 G4double *pphi, G4do 164 G4double *pphi, G4double *ptheta) const 165 { 165 { 166 G4double rand1, rand2, rand3 = 0; 166 G4double rand1, rand2, rand3 = 0; 167 G4double phi = 0; 167 G4double phi = 0; 168 G4double theta = 0; 168 G4double theta = 0; 169 G4double crossSectionValue = 0; 169 G4double crossSectionValue = 0; 170 G4double crossSectionMajorantFunctionValue = 170 G4double crossSectionMajorantFunctionValue = 0; 171 G4double maxBeta = 0; 171 G4double maxBeta = 0; 172 172 173 do { 173 do { 174 174 175 rand1 = G4UniformRand(); 175 rand1 = G4UniformRand(); 176 rand2 = G4UniformRand(); 176 rand2 = G4UniformRand(); 177 rand3 = G4UniformRand(); 177 rand3 = G4UniformRand(); 178 178 179 phi=2*pi*rand1; 179 phi=2*pi*rand1; 180 180 181 if(shellLevel == 0){ 181 if(shellLevel == 0){ 182 182 183 // Polarized Gavrila Cross-Section for K 183 // Polarized Gavrila Cross-Section for K-shell (1959) 184 theta=std::sqrt(((std::exp(rand2*std::lo 184 theta=std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta); 185 crossSectionMajorantFunctionValue = Cros 185 crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta); 186 crossSectionValue = DSigmaKshellGavrila1 186 crossSectionValue = DSigmaKshellGavrila1959(beta, theta, phi); 187 187 188 } else { 188 } else { 189 189 190 // Polarized Gavrila Cross-Section for 190 // Polarized Gavrila Cross-Section for other shells (L1-shell) (1961) 191 theta = std::sqrt(((std::exp(rand2*std:: 191 theta = std::sqrt(((std::exp(rand2*std::log(1+cBeta*pi*pi)))-1)/cBeta); 192 crossSectionMajorantFunctionValue = Cros 192 crossSectionMajorantFunctionValue = CrossSectionMajorantFunction(theta, cBeta); 193 crossSectionValue = DSigmaL1shellGavrila 193 crossSectionValue = DSigmaL1shellGavrila(beta, theta, phi); 194 194 195 } 195 } 196 196 197 maxBeta=rand3*aBeta*crossSectionMajorantFu 197 maxBeta=rand3*aBeta*crossSectionMajorantFunctionValue; 198 198 199 }while(maxBeta > crossSectionValue); 199 }while(maxBeta > crossSectionValue); 200 200 201 *pphi = phi; 201 *pphi = phi; 202 *ptheta = theta; 202 *ptheta = theta; 203 } 203 } 204 204 205 // 205 // 206 206 207 G4double G4RDPhotoElectricAngularGeneratorPola 207 G4double G4RDPhotoElectricAngularGeneratorPolarized::CrossSectionMajorantFunction(const G4double theta, const G4double cBeta) const 208 { 208 { 209 // Compute Majorant Function 209 // Compute Majorant Function 210 G4double crossSectionMajorantFunctionValue = 210 G4double crossSectionMajorantFunctionValue = 0; 211 crossSectionMajorantFunctionValue = theta/(1 211 crossSectionMajorantFunctionValue = theta/(1+cBeta*theta*theta); 212 return crossSectionMajorantFunctionValue; 212 return crossSectionMajorantFunctionValue; 213 } 213 } 214 214 215 // 215 // 216 216 217 G4double G4RDPhotoElectricAngularGeneratorPola 217 G4double G4RDPhotoElectricAngularGeneratorPolarized::DSigmaKshellGavrila1959(const G4double beta, const G4double theta, const G4double phi) const 218 { 218 { 219 219 220 //Double differential K shell cross-section 220 //Double differential K shell cross-section (Gavrila 1959) 221 221 222 G4double beta2 = beta*beta; 222 G4double beta2 = beta*beta; 223 G4double oneBeta2 = 1 - beta2; 223 G4double oneBeta2 = 1 - beta2; 224 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 224 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 225 G4double oneBeta2_to_3_2 = std::pow(oneBeta2 225 G4double oneBeta2_to_3_2 = std::pow(oneBeta2,1.5); 226 G4double cosTheta = std::cos(theta); 226 G4double cosTheta = std::cos(theta); 227 G4double sinTheta2 = std::sin(theta)*std::si 227 G4double sinTheta2 = std::sin(theta)*std::sin(theta); 228 G4double cosPhi2 = std::cos(phi)*std::cos(ph 228 G4double cosPhi2 = std::cos(phi)*std::cos(phi); 229 G4double oneBetaCosTheta = 1-beta*cosTheta; 229 G4double oneBetaCosTheta = 1-beta*cosTheta; 230 G4double dsigma = 0; 230 G4double dsigma = 0; 231 G4double firstTerm = 0; 231 G4double firstTerm = 0; 232 G4double secondTerm = 0; 232 G4double secondTerm = 0; 233 233 234 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe 234 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) * 235 (sinTheta2 * cosPhi2)/std::pow(o 235 (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)* 236 (1-sqrtOneBeta2)/(4*oneBeta2_to_ 236 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3); 237 237 238 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st 238 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) * 239 (4*beta2/sqrtOneBeta2 * sinThet 239 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2 240 - 4*(1-sqrtOneBeta2)/oneBeta2 * 240 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta 241 + 4*beta2*(1-sqrtOneBeta2)/oneB 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*one 242 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 + 243 (1-sqrtOneBeta2)/oneBeta2_to_3_ 243 (1-sqrtOneBeta2)/oneBeta2_to_3_2 * cosTheta - beta * (1-sqrtOneBeta2)/oneBeta2_to_3_2); 244 244 245 dsigma = ( firstTerm*(1-pi*fine_structure_co 245 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) ); 246 246 247 return dsigma; 247 return dsigma; 248 } 248 } 249 249 250 // 250 // 251 251 252 G4double G4RDPhotoElectricAngularGeneratorPola 252 G4double G4RDPhotoElectricAngularGeneratorPolarized::DSigmaL1shellGavrila(const G4double beta, const G4double theta, const G4double phi) const 253 { 253 { 254 254 255 //Double differential L1 shell cross-section 255 //Double differential L1 shell cross-section (Gavrila 1961) 256 256 257 G4double beta2 = beta*beta; 257 G4double beta2 = beta*beta; 258 G4double oneBeta2 = 1-beta2; 258 G4double oneBeta2 = 1-beta2; 259 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 259 G4double sqrtOneBeta2 = std::sqrt(oneBeta2); 260 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1 260 G4double oneBeta2_to_3_2=std::pow(oneBeta2,1.5); 261 G4double cosTheta = std::cos(theta); 261 G4double cosTheta = std::cos(theta); 262 G4double sinTheta2 =std::sin(theta)*std::sin 262 G4double sinTheta2 =std::sin(theta)*std::sin(theta); 263 G4double cosPhi2 = std::cos(phi)*std::cos(ph 263 G4double cosPhi2 = std::cos(phi)*std::cos(phi); 264 G4double oneBetaCosTheta = 1-beta*cosTheta; 264 G4double oneBetaCosTheta = 1-beta*cosTheta; 265 265 266 G4double dsigma = 0; 266 G4double dsigma = 0; 267 G4double firstTerm = 0; 267 G4double firstTerm = 0; 268 G4double secondTerm = 0; 268 G4double secondTerm = 0; 269 269 270 firstTerm = sinTheta2*cosPhi2/std::pow(oneBe 270 firstTerm = sinTheta2*cosPhi2/std::pow(oneBetaCosTheta,4)-(1 - sqrtOneBeta2)/(2*oneBeta2) 271 * (sinTheta2 * cosPhi2)/std::po 271 * (sinTheta2 * cosPhi2)/std::pow(oneBetaCosTheta,3) + (1-sqrtOneBeta2)* 272 (1-sqrtOneBeta2)/(4*oneBeta2_to_ 272 (1-sqrtOneBeta2)/(4*oneBeta2_to_3_2) * sinTheta2/std::pow(oneBetaCosTheta,3); 273 273 274 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(st 274 secondTerm = std::sqrt(1 - sqrtOneBeta2)/(std::pow(2.,3.5)*beta2*std::pow(oneBetaCosTheta,2.5)) * 275 (4*beta2/sqrtOneBeta2 * sinThet 275 (4*beta2/sqrtOneBeta2 * sinTheta2*cosPhi2/oneBetaCosTheta + 4*beta/oneBeta2 * cosTheta * cosPhi2 276 - 4*(1-sqrtOneBeta2)/oneBeta2 * 276 - 4*(1-sqrtOneBeta2)/oneBeta2 *(1+cosPhi2) - beta2 * (1-sqrtOneBeta2)/oneBeta2 * sinTheta2/oneBetaCosTheta 277 + 4*beta2*(1-sqrtOneBeta2)/oneB 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*one 278 + (1-sqrtOneBeta2)/(4*beta2*oneBetaCosTheta*oneBetaCosTheta) * (beta/oneBeta2 - 2/oneBeta2 * cosTheta * cosPhi2 + 279 (1-sqrtOneBeta2)/oneBeta2_to_3_ 279 (1-sqrtOneBeta2)/oneBeta2_to_3_2*cosTheta - beta*(1-sqrtOneBeta2)/oneBeta2_to_3_2); 280 280 281 dsigma = ( firstTerm*(1-pi*fine_structure_co 281 dsigma = ( firstTerm*(1-pi*fine_structure_const/beta) + secondTerm*(pi*fine_structure_const) ); 282 282 283 return dsigma; 283 return dsigma; 284 } 284 } 285 285 286 G4double G4RDPhotoElectricAngularGeneratorPola 286 G4double G4RDPhotoElectricAngularGeneratorPolarized::GetMax(const G4double arg1, const G4double arg2) const 287 { 287 { 288 if (arg1 > arg2) 288 if (arg1 > arg2) 289 return arg1; 289 return arg1; 290 else 290 else 291 return arg2; 291 return arg2; 292 } 292 } 293 293 294 // 294 // 295 295 296 G4RotationMatrix G4RDPhotoElectricAngularGener 296 G4RotationMatrix G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronRotationMatrix(const G4ThreeVector& direction, 297 const G4ThreeVector 297 const G4ThreeVector& polarization) const 298 { 298 { 299 G4double mK = direction.mag(); 299 G4double mK = direction.mag(); 300 G4double mS = polarization.mag(); 300 G4double mS = polarization.mag(); 301 G4ThreeVector polarization2 = polarization; 301 G4ThreeVector polarization2 = polarization; 302 const G4double kTolerance = 1e-6; 302 const G4double kTolerance = 1e-6; 303 303 304 if(!(polarization.isOrthogonal(direction,kTo 304 if(!(polarization.isOrthogonal(direction,kTolerance)) || mS == 0){ 305 G4ThreeVector d0 = direction.unit(); 305 G4ThreeVector d0 = direction.unit(); 306 G4ThreeVector a1 = SetPerpendicularVector( 306 G4ThreeVector a1 = SetPerpendicularVector(d0); 307 G4ThreeVector a0 = a1.unit(); 307 G4ThreeVector a0 = a1.unit(); 308 G4double rand1 = G4UniformRand(); 308 G4double rand1 = G4UniformRand(); 309 G4double angle = twopi*rand1; 309 G4double angle = twopi*rand1; 310 G4ThreeVector b0 = d0.cross(a0); 310 G4ThreeVector b0 = d0.cross(a0); 311 G4ThreeVector c; 311 G4ThreeVector c; 312 c.setX(std::cos(angle)*(a0.x())+std::sin(a 312 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x()); 313 c.setY(std::cos(angle)*(a0.y())+std::sin(a 313 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y()); 314 c.setZ(std::cos(angle)*(a0.z())+std::sin(a 314 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z()); 315 polarization2 = c.unit(); 315 polarization2 = c.unit(); 316 mS = polarization2.mag(); 316 mS = polarization2.mag(); 317 }else 317 }else 318 { 318 { 319 if ( polarization.howOrthogonal(directio 319 if ( polarization.howOrthogonal(direction) != 0) 320 { 320 { 321 polarization2 = polarization - polarizatio 321 polarization2 = polarization - polarization.dot(direction)/direction.dot(direction) * direction; 322 } 322 } 323 } 323 } 324 324 325 G4ThreeVector direction2 = direction/mK; 325 G4ThreeVector direction2 = direction/mK; 326 polarization2 = polarization2/mS; 326 polarization2 = polarization2/mS; 327 327 328 G4ThreeVector y = direction2.cross(polarizat 328 G4ThreeVector y = direction2.cross(polarization2); 329 329 330 G4RotationMatrix R(polarization2,y,direction 330 G4RotationMatrix R(polarization2,y,direction2); 331 return R; 331 return R; 332 } 332 } 333 333 334 void G4RDPhotoElectricAngularGeneratorPolarize 334 void G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronGetMajorantSurfaceAandCParameters(const G4int shellLevel, const G4double beta,G4double *majorantSurfaceParameterA, G4double *majorantSurfaceParameterC) const 335 { 335 { 336 // This member function finds for a given sh 336 // This member function finds for a given shell and beta value of the outgoing electron the correct Majorant Surface parameters 337 337 338 G4double aBeta,cBeta; 338 G4double aBeta,cBeta; 339 G4double bMin,bStep; 339 G4double bMin,bStep; 340 G4int indexMax; 340 G4int indexMax; 341 G4int level = shellLevel; 341 G4int level = shellLevel; 342 if(shellLevel > 1) level = 1; // protection 342 if(shellLevel > 1) level = 1; // protection since only K and L1 polarized double differential cross-sections were implemented 343 343 344 bMin = betaArray[0]; 344 bMin = betaArray[0]; 345 bStep = betaArray[1]; 345 bStep = betaArray[1]; 346 indexMax = (G4int)betaArray[2]; 346 indexMax = (G4int)betaArray[2]; 347 const G4double kBias = 1e-9; 347 const G4double kBias = 1e-9; 348 348 349 G4int k = (G4int)((beta-bMin+kBias)/bStep); 349 G4int k = (G4int)((beta-bMin+kBias)/bStep); 350 350 351 if(k < 0) 351 if(k < 0) 352 k = 0; 352 k = 0; 353 if(k > indexMax) 353 if(k > indexMax) 354 k = indexMax; 354 k = indexMax; 355 355 356 if(k == 0) 356 if(k == 0) 357 aBeta = GetMax(aMajorantSurfaceParameterTa 357 aBeta = GetMax(aMajorantSurfaceParameterTable[k][level],aMajorantSurfaceParameterTable[k+1][level]); 358 else if(k==indexMax) 358 else if(k==indexMax) 359 aBeta = GetMax(aMajorantSurfaceParameterTa 359 aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]); 360 else{ 360 else{ 361 aBeta = GetMax(aMajorantSurfaceParameterTa 361 aBeta = GetMax(aMajorantSurfaceParameterTable[k-1][level],aMajorantSurfaceParameterTable[k][level]); 362 aBeta = GetMax(aBeta,aMajorantSurfaceParam 362 aBeta = GetMax(aBeta,aMajorantSurfaceParameterTable[k+1][level]); 363 } 363 } 364 364 365 if(k == 0) 365 if(k == 0) 366 cBeta = GetMax(cMajorantSurfaceParameterTa 366 cBeta = GetMax(cMajorantSurfaceParameterTable[k][level],cMajorantSurfaceParameterTable[k+1][level]); 367 else if(k == indexMax) 367 else if(k == indexMax) 368 cBeta = GetMax(cMajorantSurfaceParameterTa 368 cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]); 369 else{ 369 else{ 370 cBeta = GetMax(cMajorantSurfaceParameterTa 370 cBeta = GetMax(cMajorantSurfaceParameterTable[k-1][level],cMajorantSurfaceParameterTable[k][level]); 371 cBeta = GetMax(cBeta,cMajorantSurfaceParam 371 cBeta = GetMax(cBeta,cMajorantSurfaceParameterTable[k+1][level]); 372 } 372 } 373 373 374 *majorantSurfaceParameterA = aBeta; 374 *majorantSurfaceParameterA = aBeta; 375 *majorantSurfaceParameterC = cBeta; 375 *majorantSurfaceParameterC = cBeta; 376 376 377 } 377 } 378 378 379 379 380 // 380 // 381 G4ThreeVector G4RDPhotoElectricAngularGenerato 381 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::PhotoElectronComputeFinalDirection(const G4RotationMatrix& rotation, const G4double theta, const G4double phi) const 382 { 382 { 383 383 384 //computes the photoelectron momentum unitar 384 //computes the photoelectron momentum unitary vector 385 G4double px = std::cos(phi)*std::sin(theta); 385 G4double px = std::cos(phi)*std::sin(theta); 386 G4double py = std::sin(phi)*std::sin(theta); 386 G4double py = std::sin(phi)*std::sin(theta); 387 G4double pz = std::cos(theta); 387 G4double pz = std::cos(theta); 388 388 389 G4ThreeVector samplingDirection(px,py,pz); 389 G4ThreeVector samplingDirection(px,py,pz); 390 390 391 G4ThreeVector outgoingDirection = rotation*s 391 G4ThreeVector outgoingDirection = rotation*samplingDirection; 392 return outgoingDirection; 392 return outgoingDirection; 393 } 393 } 394 394 395 // 395 // 396 396 397 void G4RDPhotoElectricAngularGeneratorPolarize 397 void G4RDPhotoElectricAngularGeneratorPolarized::PrintGeneratorInformation() const 398 { 398 { 399 G4cout << "\n" << G4endl; 399 G4cout << "\n" << G4endl; 400 G4cout << "Polarized Photoelectric Angular G 400 G4cout << "Polarized Photoelectric Angular Generator" << G4endl; 401 G4cout << "PhotoElectric Electron Angular Ge 401 G4cout << "PhotoElectric Electron Angular Generator based on the general Gavrila photoelectron angular distribution" << G4endl; 402 G4cout << "Includes polarization effects for 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-se 403 G4cout << "For higher shells the L1 cross-section is used." << G4endl; 404 G4cout << "(see Physics Reference Manual) \n 404 G4cout << "(see Physics Reference Manual) \n" << G4endl; 405 } 405 } 406 406 407 G4ThreeVector G4RDPhotoElectricAngularGenerato 407 G4ThreeVector G4RDPhotoElectricAngularGeneratorPolarized::SetPerpendicularVector(const G4ThreeVector& a) const 408 { 408 { 409 G4double dx = a.x(); 409 G4double dx = a.x(); 410 G4double dy = a.y(); 410 G4double dy = a.y(); 411 G4double dz = a.z(); 411 G4double dz = a.z(); 412 G4double x = dx < 0.0 ? -dx : dx; 412 G4double x = dx < 0.0 ? -dx : dx; 413 G4double y = dy < 0.0 ? -dy : dy; 413 G4double y = dy < 0.0 ? -dy : dy; 414 G4double z = dz < 0.0 ? -dz : dz; 414 G4double z = dz < 0.0 ? -dz : dz; 415 if (x < y) { 415 if (x < y) { 416 return x < z ? G4ThreeVector(-dy,dx,0) : G 416 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); 417 }else{ 417 }else{ 418 return y < z ? G4ThreeVector(dz,0,-dx) : G 418 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); 419 } 419 } 420 } 420 } 421 421